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In this paper we develop a gapless theory of BEC which can be applied to both trapped and 
homogeneous gases at zero and finite temperature. The starting point for the theory is the second 
quantized, many-body Hamiltonian for a system of structureless bosons with pairwise interactions. 
A number conserving approach is used to rewrite this Hamiltonian in a form which is approximately 
quadratic with higher order cubic and quartic terms. The quadratic part of the Hamiltonian can be 
diagonalized exactly by transforming to a quasiparticle basis, while requiring that the condensate 
satisfy the Gross-Pitaevskii equation. The non-quadratic terms are assumed to have a small effect 
and are dealt with using first and second order perturbation theory. The conventional treatment of 
these terms, based on factorization approximations, is shown to be inconsistent. 

Infra-red divergences can appear in individual terms of the perturbation expansion, but we show 
analytically that the total contribution beyond quadratic order is finite. The resulting excitation 
spectrum is gapless and the energy shifts are small for a dilute gas away from the critical region, 
justifying the use of perturbation theory. Ultra-violet divergences can appear if a contact potential 
is used to describe particle interactions. We show that the use of this potential as an approximation 
to the two-body T-matrix leads naturally to a high-energy renormalization. 

The theory developed in this paper is therefore well-defined at both low and high energy and 
provides a systematic description of Bose-Einstein condensation in dilute gases. It can therefore 
be used to calculate the energies and decay rates of the excitations of the system at temperatures 
approaching the phase transition. 



I. INTRODUCTION 

Dilute Bose-condensed gases provide a rare example of an interacting, many-body system for which a quantitative, 
microscopic analysis is possible at finite temperature. For this reason such systems have long been the subject of 
theoretical study 0| and much research was done in the 1950's and 60's with a view to understanding the properties of 
liquid Helium. However, the strong interactions which are present in that system mask the purely quantum statistical 
effects of condensation, and it was not until 1995 with the first production of BEC in trapped gases that a quantitative 
comparison of theory and experiment became possible Q . A large number of experiments on dilute gas BECs now 
exist and a wide range of properties have been measured with high precision. These systems therefore provide a very 
important model for testing the application of quantum field theoretical techniques to coherent, many-body systems 
at finite temperature. 

A wide variety of different (but related) approaches to the theory of BEC were developed for the homogeneous case 
and have now been extended to trapped gases. The first quantitative analysis was given by Bogoliubov in 1947 [Q and 
was based on approximating the Hamiltonian of the system by one which is quadratic and hence can be diagonalized 
exactly. The extension of this method to higher order calculations and trapped gases has become known as the 
Hartree-Fock-Bogoliubov (HFB) theory Bogoliubov obtained his quadratic Hamiltonian by using a description 
of BEC in terms of spontaneous symmetry breaking. It has recently been shown, however, that this is not necessary 
and the same results can be obtained using a number conserving approach We will follow such an approach in 

this paper. 

In 1958, Beliaev introduced a theory of BEC at T = based on the Green's functions of quantum field theory 
This approach was further developed by Popov and Fadeev ^,|l^ who applied it to calculations at finite temperature, 
and recently Fedichev and Shlyapnikov |l^] have extended the calculations to higher order and to trapped gases. A 
similar approach (but restricted to the homogeneous limit) has also been described by Shi and Griffin | |T^ . 

Beliaev's approach was also pursued by Hugcnholtz and Pines p^ , who derived the result now known as the 
Hugenholtz-Pines theorem. This theorem shows that the energy spectrum of a Bose gas is gapless, which means that 



1 



the energy of an excitation tends to zero as its momentum tends to zero. Q This is an exact result which puts a strong 
constraint on any theoretical description of BEC, and can be seen as a general consequence of spontaneous symmetry 
breaking, as shown by Goldstone . 

An alternative approach to the theory of the dilute Bose gas was introduced by Lee, Huang and Yang [|5|- jlj] and 
was based on the use of the pscudopotcntial as a means of describing low-energy collisions. Higher order calculations 
using this method were given by Wu , who calculated the ground state energy, and by Mohling and Sirlin jl^] and 
Mohling and Morita pO| , who calculated the excitation spectrum at zero and finite temperature. For some reason 
these last two papers are not commonly cited in the literature, even though they give an explicit demonstration of 
the fact that consistent calculations beyond the quadratic Hamiltonian of Bogoliubov produce finite, gapless results. 
The theory described in this paper is the logical extension of the work of Mohling, Sirlin and Morita to the case of 
an inhomogeneous gas. 

In 1959, a number conserving variational approach which included the effect of all pair correlations in a dilute 
Bose gas was introduced by Girardeau and Arnowitt [Q. However, the direct application of the variational principle 
violates the Hugenholtz-Pines theorem by predic ting a gap in the low-energy excitation spectrum (in fact this approach 
reproduces the HFB theory, as discussed in Sec. IV E ). It was later shown by Takano that the inclusion of cubic 
terms in the Hamiltonian removes this gap, a result which is consistent with the work of Mohling and Sirlin |l8| ] 
and this paper. More recently, the variational approach has been used by Bijlsma and Stoof ||2^ who considered the 
properties of homogeneous Bose gases in two and three dimensions. By combining the variational method with an 
effective Hamiltonian these authors were able to obtain a gapless theory. This approach is essentially equivalent to a 
gapless extension of the HFB theory which we will discuss in Sec. VIH. 

The low-temperature properties of trapped gases BECs are well-described by the Gross-Pitaevskii equation (GPE) 
p3| , and properties of the condensate such as its size, shape and energy can be determined from the static solutions 
to this equation p^ , p5t . The excitation energies can be found by linearization around these static solutions |2^, and 
the results of these calculations are in good agreement with experiment at low temperature |27 2^. An extension of 
the linear response approach which includes the effect of the noncondensate has recently been developed by Rusch 
and Burnett P9[ . 

Current numerical calculations of the properties of BEC at finite temperature are usually based on the Popov 
approximation to the HFB theory 1^ . The results of these calculations are in excellent agreement with experiment 
for temperatures below about 60% of the critical temperature [ pOpH |. There is a significant discrepancy at higher 
temperatures, however, and this has motivated the development of theories which go beyond HFB-Popov. In addition 
there are a number of theoretical difficulties associated with the HFB approach, such as the presence of divergences 
and the prediction of a gap in the excitation spectrum at low-energy (see below). These difficulties are avoided in 
the Popov approximation where the term responsible for these effects is neglected. This is a rather ad hoc procedure, 
however, and there is a need for a systematic approach which goes beyond the HFB theory in a consistent manner. In 
this paper we will develop such an approach and show that it leads to a theory which is gapless and free of divergences. 
It is hoped that the numerical implementation of this theory to trapped gases will yield improved agreement with 
experiment at high temperatures. 



A. Problems with HFB 



It is well-known that calculations beyond the quadratic Hamiltonian of Bogoliubov using the HFB theory encounter 
difhculties associated with the presence of infra-red and ultra-violet divergences and the appearance of a gap in the 
excitation spectrum, in violation of the Hugenholtz-Pines theorem. Ultra-violet divergences arise if atomic interactions 
are approximated by a contact potential via 

V{r)^UoS{r), C/q = , (1) 

m 

where a is the s-wave scattering length and m is the atomic mass. The contact potential is the lowest order approxi- 
mation to the pseudopotential of Huang and Yang flB] , and the problem of ultra-violet divergences can be avoided by 



^An exception is the charged Bose gas, for which the Hugenholtz-Pines theorem is satisfied but nonetheless there is a gap 
in the excitation spectrum at low energy because the Fourier transform of the Coulomb potential is proportional to as 
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using a better approximation,g as in Ref. [Q. This has the disadvantage, however, that it produces a Hamiltonian 
which is not Hermitian. It is also not sufficiently general for our purposes as it does not allow the theory to be used 
for the study of two-dimensional or charged gases. We will therefore develop the theory using V{r) and the results 
we will obtain can be applied to a wide range of physical systems. For the special case of neutral atoms in three 
dimensions, we will show how the contact potential can be introduced in a manner which does not lead to ultra-violet 
divergences. 

Infra-red divergences also appear in the HFB theory. An example is provided by the perturbative shifts to the 
quasiparticle energies which scale with momentum fc as ~ 1/fc as fc — > (see Sec. VI). One way of dealing with 
these difficulties is to use a description of the system in terms of density and phase variables, as shown by Popov Q . 
We will show instead that the infra-red divergences are spurious, in the sense that a consistent treatment of higher 
order terms leads to finite predictions for physical quantities, and also to a gapless excitation spectrum. 



B. Outline of paper 

The development of the theory starts in Sec. |^ from the second quantized, many-body Hamiltonian for a system of 
identical, structureless bosons. Interactions between the atoms are assumed to correspond to binary collisions and are 
described by a general interatomic potential V{r). When a condensate is present, the Hamiltonian can be rewritten 
in the form of Eqs. ([^)-(|l^) and is quadratic to leading order with higher order terms which are cubic and quartic 
in noncondensate operators. This BEC form of the Hamiltonian is usually obtained using spontaneous symmetry 
breaking, but instead we will use arguments based on the conservation of particle number following the approach of 
Refs. [g-g. 

The BEC Hamiltonian allows a systematic calculation of the properties of dilute gas BECs, because any quadratic 
Hamiltonian can be diagonalized exactly, while the non-quadratic terms should be small and can be dealt with 



perturbatively. The exact diagonalization of the quadratic part is the subject of Sec. [II, and is achieved using 
a transformation to a quasiparticle basis, with the condensate wave function determined by the time-independent 
GPE. The quasiparticle energies and transformation coefhcients are calculated from the Bogoliubov-de Gennes (BdG) 
equations. 

The cubic and quartic terms in the Hamiltonian produce shifts and widths of the quasiparticle energies and become 
progressively more important as the temperature increases and a significant noncondensate fraction appears. These 
terms are treated using first and second order perturbation theory in a quasiparticle basis. Sec. ^ is concerned 
with the explicit calculation of the perturbative expressions in terms of the quasiparticle energies and transformation 
coefficients. We also show that the conventional treatment of the non-quadratic Hamiltonian in terms of factorization 
approximations is inconsistent. 

In Sec. ^ we discuss the physical interpretation of the higher order terms obtained from perturbation theory. We 
show that they modify the description of particle scattering and upgrade the bare interatomic potential V{r) to a 
T-matrix. This allows us to rewrite all interaction matrix elements in terms of the two-body T-matrix which describes 
particle collisions in a vacuum. It is this T-matrix, not V{r), which can be approximated by a contact potential for 
low-energy scattering. The difference between the two-body T-matrix and V'(r) naturally provides the renormalization 
which is required to remove ultra-violet divergences in the theory. 

The problem of infra-red divergences in the theory of BEC is a feature of the homogeneous limit since in a trapped 



gas there is a natural low-energy cut-off determined by the size of the trap. In Sec. VI we therefore specialize to this 
limit and present explicit results for the energies and lifetimes of the quasiparticle excitations. We show analytically 
that, although infra-red divergences do occur in individual terms of the perturbative expressions, the total contribution 
beyond quadratic order is finite. In addition, the excitation spectrum is gapless as required by the Hugenholtz-Pines 
theorem. These results are valid at both zero and finite temperature and prove that the theory developed in this 
paper provides a consistent description of BEC beyond the approximation of a quadratic Hamiltonian. 



In Sec. VII we discuss the expected range of validity of the perturbative treatment of non-quadratic terms. We 
show that this approach is justified for a dilute gas, but that in the homogeneous limit it must fail at the critical 
point. For a trapped gas, however, it is possible that the theory will remain valid even in the region of the phase 



transition. We conclude in Sec. VIII with a description of a gapless extension to the usual HFB theory |M] which can 



■^Specifically the correct form of the pseudopotential is V{r) — Uo 5{r){d/dr)r + 0[{ka)^], where r = ri — r2 is the separation 
of the two interacting atoms and k is their relative momentum. 
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be used as an approximation to the full theory of this paper. This approach has recently been shown to give improved 
agreement with experiment at temperatures approaching the critical point . 

II. THE BEC HAMILTONIAN 

The many-body Hamiltonian for a system of structureless bosons with pairwise interactions can be written in the 
usual second quantized formalism as 

H = J2 H^^^3 + \Y. mV'\km)a\a]akara. (2) 

ij ijkm 

The operators a| and respectively create or annihilate a particle from the basis state with wave function Q (r) and 
obey the usual Bose commutation relations. The matrix elements H^j' are defined by 

h:^^ Id^^cw^^^ow, (3) 

where the single-particle Hamiltonian H'^p is 

H''-^^^' + VTr.,{r), (4) 

and VTrap(r) is the magnetic potential that confines the atoms. For current experimental configurations this can be 
taken to be harmonic, although the formalism can also be used to describe a homogeneous gas if the trap potential 
is set to zero. We will denote the eigenvalues of H^p by e^^ and the corresponding eigenvectors (normalized to 1) by 
(^^^{r). In general, these functions do not correspond to the basis functions used in Eq. (^) [i.e. Ci(r) 7^ Cf^(r)] and 
consequently H-j is not diagonal (see below). 

The interaction matrix elements in Eq. (|^) are symmetrized and are defined by 

{ij\V'\km) = i \{ij\V\km) + {ji\V\km)] , (5) 
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where 



{ij\V\km) = J d^n d'r^ {C(ri)C;(r2)^(r)a(r2)Cm(ri)} , (r = ri - ra), (6) 

and V{r) is the bare interatomic potential. 

The Hamiltonian of Eq. (^) is written in the basis provided by the orthonormal functions {Ci(r)}. One of these, 
Co (r) , describes the condensate and we will show in Sec. HI that it is a solution of the time- independent GPE [E3| . 



The remaining functions Ci(r) (i 7^ 0) form a complete set orthogonal to the condensate. The choice of these functions 



is a matter of convenience, and in Sec. HI we will transform to the quasiparticle amplitudes Ui(r) and Vi(r) which are 
determined by the Bogoliubov-de Gennes (BdG) equations. Two different choices for the basis functions will lead to 
the same quasiparticle amplitudes, but via a different transformation. We note, however, that for a trapped gas, the 
functions Q (r) can not be chosen to be eigenstates of the single-particle Hamiltonian because they must be orthogonal 
to the condensate and this is not an eigenstate of H'^p in a trapQ This issue does not arise in the homogeneous limit 



where the functions (i{r) can be taken to be the usual plane waves (see Sec. VI). 

When a condensate is present the Hamiltonian of Eq. (^ can be rewritten in a form which allows a systematic 
approximation scheme to be developed. This can be achieved using a number conserving approach, following the 
arguments of Girardeau and Arnowitt , Gardiner ||] and Castin and Dum . The details of these arguments are 
given in Ref. and we will merely quote the results here. 

Number conservation is achieved by defining the pair operators cti = P^di, where /3o = (-^0 + l)~^^^flo and Nq 
is the number operator for the condensate. The {a^} clearly conserve particle number, and they also satisfy Bose 



Possible choices for these functions are given in Eqs. (|l|) and 
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commutation relations exactly provided only that the subspace of total condensate depletion is ignored. For all 
practical purposes they can therefore be treated as ordinary Bose operators and in what follows we will simply denote 
them by Oj rather than ai. 

It is straightforward to show that the Hamiltonian can be rewritten in terms of the pair operators in the form |3q] 



where 



Ha 



(7) 



(8) 



h.c 



(9) 



-A%+2iVo(0^|F1J0) 



E 



^{^J\V'^\00)ala] + h.c. 



(10) 



E 



(11) 



H4 = 



k'm)aia)akam- 



(12) 



Here the average (. . .) stands for a quantum expectation value in a pure energy eigenstate (which for the purposes of this 

paper means a quasiparticle number state), while N^x is the number operator for the noncondcnsate, N^x = Si^^o '^I'^i- 

The parameter iVo is the mean condensate population and is defined by A'o = — {Nex)- In Eq- (0) we have also 
defined the parameter A by 



A 



-"00 



JVo(00|F^|00). 



(13) 



We will show in Sec. Ill that this quantity is in fact the condensate eigenvalue as calculated from the GPE. The 
number conserving approach automatically leads to this parameter entering the quadratic Hamiltonian, and we note 
that it is distinct from the related parameter fx (the chemical potential) which appears in H2 in a broken symmetry 
approach (see below). 

The number conserving Hamiltonian given above has essentially the same form as that obtained using spontaneous 
symmetry breaking. In this approach, the condensate annihilation operator oq is replaced with the number \/Nq, on 
the grounds that the commutation relations of oq can be neglected when Nq is macroscopic Since this procedure 
docs not conserve particle number, the grand canonical Hamiltonian H — iiN is diagonalized rather than simply H, 
the chemical potential /i being chosen (as a function of temperature) so that the mean number of particles is constant. 

The principle differences between the number conserving and broken symmetry approaches are in the interpretation 
of the operators, the replacement of the condensate eigenvalue A with the chemical potential in the quadratic 
Hamiltonian and the fact that the condensate population Nq is not an independent parameter in a number conserving 
approach, but must be determined self-consistently by iVo — N ~ {Nex)- In addition, the BEG Hamiltonian of 
Eqs. (0)-(|2| describes the energy of the system rather than the 'free energy' H — ^N which is obtained in spontaneous 
symmetry breaking. The relation between the two approaches is discussed in more detail in Ref . p5| . 

The Hamiltonian of Eqs. (^-(|2|) is the basis of a systematic treatment of a Bose condensed system. The presence 
of a condensate means that Nq is a large number and we can therefore group terms according to the powers of the 
condensate population they contain. The dominant terms should correspond to the quadratic Hamiltonian Hq — 
Hq + Hi + H2 which can be diagonalized exactly, while the non-quadratic terms should have a small effect and can be 
dealt with perturbatively. The results of this proce dure are discussed in the following sections, while a discussion of 
the validity of perturbation theory is deferred to Sec. VII when the theory has been developed further. It will turn out. 
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however, that at T = the use of perturbation theory is justified if na^ ^ 1, where n = N/Q. is the particle number 
density. This is the usual dilute gas criterion and is well-satisfied in current experiments for which na^ ~ lO"'^. At 
finite temperature, this condition is replaced by the requirement (fcf,T/nn?7n).(nna^)"^^^ ^ 1, where ng = Nq/Q. is the 
condensate density. This is the same result as that obtained in Ref. 



III. THE QUADRATIC HAMILTONIAN 



In this section we discuss the exact diagonalization of the quadratic Hamiltonian. The first stage is described in 
Sec. [II A and involves the minimization of the energy functional iJp. This requires the condensate wave function to be 
a solution of the GPE, which in turn means that the linear Hamiltonian Hi is identically zero. The diagonalization of 
i?2 is described in Sec. [II B , and is achieved by transforming to a quasiparticle basis. The quasiparticle energies and 
wave functions are determined by the Bogoliubov-de Gennes equations. These equations are rewritten in the position 
representation in Sec. [II C, where we also discuss the issue of orthogonality of the excitations to the cond ensate. 
Although the development of the theory in this section applies to pure states, we conclude in Sec. [II D with a 
discussion of how thermal averages can be obtained. 



A. The GPE 

The first stage in the diagonalization of the quadratic Hamiltonian is the minimization of ffoi which is a functional of 
the condensate wave function Co(r). The eigenstate solutions for the condensate can therefore be found by functional 
differentiation of this Hamiltonian with respect to Q{v). This leads immediately to the time- independent GPE which 
determines the condensate wave function ||23| . In basis space notation this can be written as the set of equations (for 
all A;) 



HZ + No{k{)\V'\m) = \5ko- 



(14) 



The GPE is usually written in the position representation using the contact potential approximation of Eq. (|l|) 
this case Eq. (Oh takes the more familiar form 



In 



^V\o(r) 
Im 



^Trap(r)Co(r) +iVo(7o|Co(r)PCo(r) = ACo(r). 



(15) 



The quantity A in Eqs. (|lj) and ( [Tsl ) is introduced as the Lagrange multiplier for the constraint on the normalization 
(/ d^T |Co(r)P = 1) when we perform the functional differentiation. Comparison with Eq. ( p^ shows that it is the same 
parameter as we defined earlier and which appears naturally in the quadratic Hamiltonian of Eq. ( p^ in a number 
conserving approach. Equation ( p^ clearly shows that A corresponds to the energy of a particle in the condensate. 

The basis wave functions C,i{v) which describe the noncondensate must be chosen to be orthogonal to Co(r). A 
convenient choice for these functions for a trapped gas is therefore provided by the solutions of the linear Schrodinger 
equation with an effective potential produced by the trap and the condensate density |30 



FTrap(r)+Afot/o|Co(r)|' 



'O(r), 



(16) 



where ef denotes the energy of the basis state. Comparison with Eq. ([l5|) shows that the solutions of this equation 
are indeed orthogonal to the condensate. 



B. Diagonalization of H2 



If the GPE of Eq. ( |l4|) is satisfied, then the linear Hamiltonian Hi of Eq. (g) vanishes. This is what we would 
expect since the condensate gives a functional minimum of Hq and linear variations vanish at extrema. The next 
order contributions to the energy of the system therefore come from the quadratic Hamiltonian H2 of Eq. (|o|). This 
can be written as 



Ho 



E 



(17) 
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where we have defined the quantities 



A, = H^P ~ A% + 2No{0t\V^\j0), 



No{ij\v'm. 



(18) 



We note that the contribution from A in Cij ensures that the excitation (quasiparticle) energies are measured relative 
to the condensate, while the expression X{Nex) in H2 is simply a number which shifts the zero of energy so that the 
total Hamiltonian gives the absolute energy of the noncondensate. 

The interaction terms in Eq. (^8|) represent the most important scattering processes which occur in a condensed 
gas. The term 2No{0i\V'^\j0) in £y describes direct (Hartree) and exchange (Fock) collisions of excited atoms off 
the condensate, while Aiij describes the effect of pair formation out of the condensate. The diagonalization of the 
quadratic Hamiltonian corresponds to summing over these direct, exchange and pair excitation processes to all orders 
in the interaction potential V{r) and is therefore equivalent to the Random Phase Approximation ]3^ . 

H2 can be diagonalized by transforming to a quasiparticle basis in which new operators Pi are defined by 



(19) 



We note that the quasiparticle transformation is defined in the subspace orthogonal to the condensate. If we have 
m basis states spanning this space, then the coefficients Uij and Vij form two m x m matrices U and V. The 
transformation of Eq. (|l^) is required to be canonical, which means that it preserves the commutation relations and 
leads to bosonic quasiparticles. This is the case if the matrices U and V satisfy the following orthogonality and 
symmetry conditions 



UW - VV^ = 1, 



(20) 



These conditions imply that the inverse quasiparticle transformation has the form 



(21) 



The necessary and sufficient conditions for the quasiparticle transformation to diagonalize H2 are given by the 
Bogoliubov-de Gennes (BdG) equations. These can be written as the matrix eigenvalue problem 



C M 

-M* -L* 



(22) 



where C and M are the matrices with elements and Aii 



by Up = ( U^ 



pi I 



p2, 



pm 



Y and Vp = {Vpi, Vj, 



p2, 



respectively and we have defined the vectors Up and Vp 
)'^. For obvious reasons, we will refer to the matrices 



£ and A4 as the diagonal and off-diagonal elements of the BdG equations respectively. If Eq. 
the Hamiltonian takes the form 



2[i is satisfied, then 



h2=j: 



1 



(23) 



while the quasiparticle energies Cp are given by 



( 



C M 

-M* -C* 



(24) 



The detailed properties of the BdG equations are discussed in Ref. [Q. We note here, however, that the solutions 
come in pairs with positive and negative energies, because if there is one solution with gp > then there is also 
another with Cp' = — Cp, Upi — and Vpi — u*. The problem of a 'missing eigenvector' may therefore appear, because 
if ep = only a single eigenvector is obtained rather than two independent ones. The solution with zero energy is 
proportional to the condensate wave function, however, so this issue does not arise in our formalism in which the BdG 
equations are written in the space orthogonal to the condensate. The solutions of Eq. (|2^) span this space and there 
is no solution with zero energy and no difficu lty w ith a missing eigenvector. A discussion of the orthogonality of the 



excitations to the condensate is given in Sec. [II C 



Despite the fact that there is no solution of the BdG equations with exactly zero energy, we will nonetheless prove 
that the excitation spectrum is gapless by showing in the homogeneous limit that the energy of a quasiparticle tends 
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to zero as its momentum tends to zero (see Sec. VI). Thus, although there is no solution to Eq. ( p2| ) with = 0, 
there may be solutions with an arbitrarily small energy. 

The BdG equations can be solved straightforwardly in the high-energy limit (ep — > oo) where the condensate has a 
small effect. In this limit the basis wave functions can be taken to be eigenstates of H'^p [i.e. Ci(r) = Ci'^(r)], while 
the quasiparticle transformation reduces to Uij — *■ Sij and Vij —^ 0. More precisely, the leading order contributions to 



Uij, Vij and €p are 



M* 

C^..— V,,^ ( 1 ^ ' (25) 



Cpp = e;P-X + 2No{Op\V'\pO). (26) 



The result of Eq. (^6|) gives the quasiparticle energy in the Hartree-Fock approximation. 
For the further development of the theory, it is convenient to define the quantities 



Pij = {a]ai), = {cLjCLi), (27) 

with the properties pij — p*,^ and = Hji. We will refer to pij as the one-body density matrix and to as the 
anomalous average. As before (. . .) denotes an ordinary quantum expectation value in a quasiparticle number state. 
The population of the noncondensate is related to pij by 

Using Eq. ( pT| ) , pij and kij can be written in terms of the quasiparticle populations Up and transformation coefficients 
Uij and Vij as 

p/0 



p#0 



[Up^Vp) + Up,Vp*^np + Up,Vp\. (30) 



We will also be interested in the change in pij and Kij when a quasiparticle is added to some particular mode p (i.e. 
Up Up + 1), which we will denote by Apij(p) and AKij{p). These quantities are given by the expressions in square 
brackets in Eqs. (||) and (|^) 

Ap,,{p) = Up,U;^ + V;,Vp,, An,,{p) = Up,Vp) + Up,Vp\. (31) 

Eqs. (IH) and ( p9| ) show that the zero-temperature depletion of the condensate depends only on the matrix V and 
hence arises from the off-diagonal part of the BdG equations which describes pair excitation from the condensate. 

Since H2 is diagonal in a quasiparticle basis, we can find its contribution to the total energy of the system by taking 
its expectation value in a quasiparticle number state. Using Eqs. (|l^), ([l8| ) and (|2^) this gives 



E2{ni\ = {H2) =y,<\Htf + 2NomV^m p,, + ^^{^J\V^\00)K* +C.C }, (32) 



2 



where c.c is the complex conjugate and pji and Kji are calculated from Eqs. ( p9|) and ( ^ ) for some quasiparticle 
distribution {ni}. We note that A does not appear in this equation as the two contributions in Eq. (|l^) cancel when 
we take the average. The total energy of the system is given to quadratic order by 

E,{n,} = Ho{n,} + E2{n,}, (33) 



*At high energy we can replace -I- ej with e^^ + in the energy denominator of Vij, but we have written the result in this 



form for the benefit of later comparison. 
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where _ffo depends on {n^} via its dependence on Nq = N — (N^x) = N — X^i/o P"- 

The excitation energies correspond to the change in the energy of the system when a single quasiparticlc is created, 
and can therefore be calculated from the expression 



= El [{rii}, Up + l] - E., [{nj, n. 
Linearizing this equation and using Eqs. (0) and (p2[), we obtain 



E 



(34) 



(35) 



with Cij and Mij as in Eq. ( p^ ) and Apij{p) and AKij{p) given by Eq. (pT|). This result is of course the same 
expression for the energy of a quasiparticlc as can be obtained directly from the BdG equations [c.f. Eq. (p4|)]. The 
method of derivation given here is more useful, however, when we come to consider the contribution of non-quadratic 
terms. 

We note that the parameter A (contained within Cij), which is absent in the expression for E2 of Eq. (^), is 
introduced into the quasiparticlc energy as the contribution from the change in Hq when A'o changes. This term 
therefore takes into account the fact that the creation of an excitation requires the removal of particles from the 
condensate. As a result the quasiparticlc energies are measured relative to the condensate. 



C. The Position Representation 



The BdG equations are often quoted in the position representation using the contact potential and we will therefore 
rewrite some of the previous equations in this more familiar form. In so doing we will address the issue of the 
orthogonality of the excitations to the condensate which arises in this representation. 

The spatial representation of the quasiparticlc transformation coefficients Uij and Vij of Eq. ( pj| ) is given by the 
functions Ui(r) and Vi{r), defined by 

u,(r) = J2 UM^). v*{v) = ^^K^- W- (36) 

The orthonormality and symmetry conditions of Eq. (|20|), which ensure that the quasiparticlc transformation is 
canonical, become the integral relations 

' ^ (37) 
d^r lui{r)vj{r) ~ Uj{r)vi{r)'> = 0. 

Using the contact potential, the BdG equations can be written as 

C{r)uj{r)+M{r)vj{r) = ejUj(r) + CjCo(r), 
£(r)«,(r)+M*(r)u,(r) = -e,z;,(r)+c,Co(r), 

where 

C{r)^H'P-X + 2NoUoMr)\\ M{r) ^ NoUoC^{r). (39) 

and 

'd^r [/o7Vo|Co(r)P \Q{r)u,{r) + Co(r)«,(r)l . (40) 



The parameters cj ensure that the solutions to Eq. ( pSf) with ej 7^ are orthogonal to the condensate . We note 
that the equations also have a solution at zero energy with Uj{r) = Co(r) and Vj{r) = — CQ(r) because the theory is 
gapless. This does not correspond to a solution of Eq. (^), however, because it is not orthogonal to the condensate. 

The BdG equations are usually written in the position representation without the parameters cj on the right hand 
side, i.e. as 
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C{r)uj(r)+M{r)ij{r) = ejUj(r), 

(41) 



These equations have the same eigenvalue solutions as Eq. ( p8D and they also have a solution with zero energy of the 
same form as that given above. Substituting this into the orthogonality and symmetry conditions of Eq. ( |37| ) shows that 

the solutions for ej ^ are orthogonal to the condensate in the generahzed sense / £'v Q{v)Uj(y) + (o(r)wj(r) = 0. 

The integrals / d'^r Cq (r)t2j(r) and J d'^r Co(r)wj(r) are not separately zero, however, and so the quasiparticle functions 
Uj(y) and ^^(r) are not individually orthogonal to the condensate. Consequently, they do not correspond to the 
functions Uj(r) and u_, (r) defined in Eq. (^6h. These functions can be obtained from the solutions of Eq. (^Tj) by 
removing the projection onto the condensatm 

= - CjCo(r)/ej, Vj(v) = f}j(r) + c^CoW/ej, (42) 

for tj 7^ 0. The significance of this result numerically is that it allows the BdG equations to be solved using a basis 
set which consists of the eigenstates of (which are not orthogonal to the condensate). 

We will complete this section by giving the (local) spatial representations of the one-body density matrix pij and 
the anomalous average which were defined in Eq. (pTj). Denoting these by Pea;(r) and respectively, we have 

Pexiy:) = ^ Cj WOWpij, K(r) = ^ GWOW'ty- (43) 

Using Eqs. (p9|), ( pO[ ) and (^) these can be written in terms of the quasiparticle amplitudes as 

p,.(r) = ^(|«p(r)p + \v,{v)\^)n, + ^^(r)^, (44) 

«(r) = ^up(r)t;;(r)(2n, + l). (45) 
The change in these quantities when np ^ + 1 will be denoted by App(r) and AKp(r). 



D. Thermal Averages 



All the averages (...) which have appeared to this point have referred to ordinary quantum expectation values 
in pure quasiparticle number states. For comparison with experiment, however, we are more interested in thermal 
averages, which should strictly be calculated via the canonical partition function Zc = X^jn;} e"''-^'^"'^ with Ei{ni\ 
as in Eq. (|33|). To a good approximation, however, they can be obtained simply by replacing the condensate and 
quasiparticle populations with their thermal averages, in which case the BdG equations and quasiparticle energies 
become dependent on temperature (via their dependence on Nq). 

The quasiparticle populations can be taken to have the usual Bose-Einstein distribution 

= ,-1,1 _ 1 (46) 

where the fugacity is defined by^z — e^^'^^^^ and is calculated from the condition that the thermal population satisfies 
{Nex) ~ N ~ No{T). The condensate population No{T) should strictly be determined from the requirement that the 
system minimizes its free energy — —ksTlogZc- The fugacity only differs significantly from unity, however, when 
the condensate population is of order 1, and so it should be sufficient simply to use the noninteracting gas result 

No{T)^z/{l~z). 



^This procedure also works for the generalized BdG equations which are used in the HFB-Popov and gapless HFB theories. 
The method fails for the full HFB theory, however, because this is not gapless (and so the condensate does not provide a zero 
energy solution of the equations). 

®The fugacity is defined here by z = e'^'''"^' rather than hy z = e'^'^ because the quasiparticle energies are measured relative 
to the condensate. 
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We note that the Bose- Einstein distribution makes use of the fact that the condensate eigenvalue and the chemical 
potential are not identical. In the limit that Nq ^ oo the difference disappears and Eq. (|4q ) reduces to the Planck 
distribution {z = 1). The distinction is important for finite systems, however, especially near the critical point where 
the Planck distribution can lead to the number of excited particles being greater than TV. 



IV. BEYOND THE QUADRATIC HAMILTONIAN 



In this section we will describe how the properties of a Bose condensed gas can be calculated beyond the approxi- 
mation of the quadratic Hamiltonian. The non-quadratic terms are expected to be small and we will therefore treat 
them using first and second order perturbation theory. This allows a calculation of the quasiparticle shifts and widths 
at zero and finite temperature. In addition, the GPE is upgraded to include the effect of the noncondensate on the 
energy and shape of the condensate. 



A. The Generalized GPE 

The GPE of Eq. ( p^ ) was obtained by requiring that the condensate wave function minimize the energy functional 
-ffo- This is a special case of the more general method of this section and is sufficient if we are only interested in the 
quadratic Hamiltonian. For calculations at higher order, however, we should include the effect of the noncondensate 
on the condensate atoms. This can be achieved by incorporating the quadratic Hamiltonian in the energy functional 
for the condensate and leads to a generalized GPE. 

The condensate wave function is determined in general by the requirement that the system minimizes its free energy 
J- — — fc^TlogZc, where Zc is the canonical partition function defined above. From the definition of J- and Zc we 
have 

where Ei is the energy of the system for some particular quasiparticle distribution and ((. . .)) denotes a thermal 
average. To obtain the GPE of Eq. (^4|) we calculated Ei using only Hq, but at this order the more accurate result 
of Eq. (|3|) is required. Using the expression for E2 from Eq. (^2|) we obtain the generalized GPE 

H^P + 7Vo(fcO|\/^|00) + J2 [2{k^\V'mPJ^ + (fcO|F^|zj>,,] = Ag4o, (48) 

where pij and Kij are to be interpreted here as thermal averages. We have denoted the generalized condensate 
eigenvalue by Xq to distinguish it from the parameter A which was introduced in Eq. (^3|). Since Xq must be real, we 
have 

J2 (00|ylzj>,, = J2 fo1^^lOO)4. (49) 

If we assume a contact potential and use the fact that the condensate wave function can be chosen to be real, then 
this result shows that the anomalous average K(r) of Eq. (|4^) can also be chosen to be real. 

If we rewrite Eq. (^ ) in the position representation using the contact potential approximation, then it takes the 
more familiar form 

- ^V'Co(r) + yT,-ap(r)Co(r) + iVoC/o|Co(r)PCo(r) + 2t/oPex(r)Co(r) + C/o«(r)Co (r) = Ac Co(r). (50) 
For a trapped gas, a convenient choice for the basis states ^^(r) is therefore provided by the solutions to the equation^ 



''To show that this gives wave functions orthogonal to the condensate we need to use the fact that Co(r), ^(r) and Q{r) can 
be chosen to be real. 



11 



2m 



V\.(r) 



VTrap(r) + iVoC/o|Co(r)P + 2C/oPe.(r)G(r) + ;7oA^(r)C(r) = ef C.(r), (51) 



which replaces Eq. (Jig). 

The use of the generahzed GPE leads to a change in the energy of the condensate from A to Xq and also to a change 
in its shape (this effect is absent in the homogeneous limit). ^ These changes in the condensate affec t th e exci ted 
atoms and produce shifts in the quasiparticle energies. These effects are discussed further in Sees. IV C and IV D . 

As a result of the additional terms in Eq. ( ^8|) relative to the ordinary GPE, the linear Hamiltonian Hi of Eq. (H) 
no longer vanishes. The remainder (which we will call AHi) is 



AHi = -VlVo J2 I ['^{k^\V'\JO)pJ^ + {kO\V'\ij)n,, 



at + h.c. }. 



(52) 



AHi has the same characteristic si ze as (^ V^o) ^-nd it can therefore be incorporated naturally into a redefinition 
of the cubic Hamiltonian (see Sec. IV D). 



B. The Cubic and Quartic Hamiltonians 



The conventional treatment of the non-quadratic terms for trapped gases is based on factorization approximations 
which reduce them to linear or quadratic forms respectively j^. These approximations are obtained by pairing 
operators in all possible ways and then replacing these pairs by their expectation values in a quasiparticle state. Thus 
we have for example 

alctjCLk 7^ PjiUk + Pkiaj + KkjO,], (53) 
+ Pmja\ak + Kmfea|a] + Kl^akCim (54) 

~ {PkiPmj + PmiPkj + ^^ijl^mk), 



where we have used the symbol -/^ to indicate that we will not be using these approximations. Equations (53) and 
( p^ ) lead to a modified quadratic Hamiltonian which can be diagonalized exactly by a (self-consistent) quasiparticle 
transformation [with pij and calculated from Eqs. ( p9| ) and (|30|)]. 
Equation (BJ) is usually justified using Wick's theorem which gives 



a|aUfca,„)) = pkiPmj + PkjPmi + K*jKmk, (55) 



while Eq. (53) is justified by analogy. In addition, these approximations lead to the same quadratic Hamiltonian as 
a variational approach ^T\ . This is not entirely satisfactory, however, and instead we will simply apply perturbation 
theory to and on the assumption that they are small compared to the quadratic Hamiltonian. This assumption 
is based on the fact that these terms contain smaller powers of the condensate population, and its validity will be 



confirmed by an explicit calculation of the energy shifts in the homogeneous limit (see Sec. VI). 

Since the condensate wave function now obeys the generalized GPE, the perturbing Hamiltonian includes the 
contribution from AHi and is given by 

Hpcrt = H3 + AHi + Hi. (56) 

We will treat this Hamiltonian using first and second order perturbation theory in a quasiparticle basis. Only H4 gives 
a non-zero contribution in first order perturbation theory, while at second order we only need to consider H3 + AHi 
by virtue of the extra factor of \/T% which this contains relative to H^i [c.f. Eqs. (Ill]), (|l2|) and (p2|)]. For the same 



We have not introduced an additional notation to describe the change in shape, so in the rest of this section the index 
which appears in matrix elements is to be interpreted as referring to a solution of the generalized GPE of Eq. ( jisl ) rather than 
the ordinary GPE of Eq. (0). 
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reason, the first order calculation with H4 is of the same order of magnitude as the second order calculation with 
+ AHi, so both must be taken into account if the calculation is to be consistent. 
In the following sections we will give the formal expressions resulting from these perturbative calculations. These 
expressions can be evaluated analytically in the homogeneous limit (see Sec. but otherwise can be incorporated 
into existing numerical codes |3^. We also show that if the first order perturbation theory on H4 is made self- 
consistent, then it leads to the same results as Eq. (jsj). The factorization approximation on H4 is therefore justified 
as long as a first order perturbative treatment is valid. However, the factorization approximation on the product of 
three operators is not justified and amounts to neglecting a number of terms which are obtained from second order 
perturbation theory. This leads to a gap in the excitation spectrum and a failure of infra-red divergences to cancel in 
physical quantities. In contrast, the perturbative treatment of the cubic Hamiltonian leads to a finite, gapless theory. 



C. First Order Perturbation Theory 



In first order perturbation theory, the energy shift to a quasiparticle number state \s) = \ni, n2, 

Epci-tis, 1) = (s|7?Port|s), 



IS 



(57) 



where the perturbing Hamiltonian Hpcrt is given in Eq. ( pq ) . Only the fourth order Hamiltonian of Ec^. (|lj) contributes 
in this expression and so we will denote this first order energy shift as E4. 

To evaluate E4 we need to calculate the expectation value of the product of four operators which appears in Eq. (|lz 
This is straightforward and the result is 



O-^djClkO'm) — PkiPmj "I" PmiPkj "I" l^ijl^mki 



(58) 



where we have left out terms which are negligible in the thermodynamic limit (the full expression is given in Ref. p5[ ) . 
This result is consistent with Wick's theorem [c.f. Eq. (p5|)], although in that casepij and are defined in terms of 
thermal averages whereas Eq. (58) has been derived for a pure state. Using Eq. (|5^) we find that E4 is given by 



E4 = {H4) = \ {W\km)[, 



PkiPmj H" PmiPkj H" ^ij ^mk\ 



(59) 



ijhm^O 



Equation ( p9| ) gives the contribution to the total energy of the system from first order perturbation theory on the 
non-quadratic Hamiltonian. For comparison with experiment, however, we are more interested in calculating the 
frequencies at which the system responds to external perturbations. These correspond to the quasiparticle energies, 
i.e. the energy required to create a quasiparticle. We are therefore interested in the change in E4 as Up Up + 1, 
which we denote by Ai?4(p) 



Ai?4(p) = i?4(ni, 7l2 . . . Tip -f 1 . . .) — i?4(ni, 712 ■ ■ - Up . . .). 



(60) 



The dependence of E4 on the quasiparticle populations is contained in pij and so we can obtain an expression for 



Aii'4 (v) simply by writing 
Eq. (ml). This gives 



Pi J + Apijip) and Ki 



AKij{p), where Apij{p) and AKijij)) are defined in 



AE4{p)^ J2 

ijkm,^0 



2{ki\V''\jm)p^kApji{p) + '^{ij\V''\k'm)KkraAK*^{p) + '^{km\V''\ij)Kl^AKj,{p) 



In the position representation (using the contact potential approximation) this becomes 

AE4{p) ^Uo J (fr [2pex(r)App(r) + K(r)AKp(r)] , 



(61) 



(62) 



where we have used the fact that k{y) can be chosen to be real. 

We must also take into account, however, that the total number of particles is fixed, so a change in the quasiparticle 
distribution must be accompanied by a change in the condensate population. This means that the creation of a 
quasiparticle via Up — s- rip -I- 1 also leads to a change in the energy of the quadratic Hamiltonian which is 
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A(Fo + ^2) = ] i^oo + NomV^m + E [m\V'\JO)pJ^ + mV^\iJ)^^,^] \ANo 

E XGS^,Ap,,{p), (63) 



where we have used Eqs. ( |4S| ) and (49) and the fact that Nq = N — (Nex)- However, the contribution to Ag from the 
ordinary GPE is included in the BdG equations of Eq. (^) and is diagonahzed exactly. The perturbative shift in the 
quasiparticle energy as a result of Eq. ( |63| ) is therefore 

AExip) = - ^ [Ag - X]S^jApj,ip) = -[Ag - X]A{Ne.){p), (64) 
where A{Nex) (p) is the change in the population of the noncondensate when a quasiparticle is created in mode p. 



We note that for a trapped gas Ag — A is not simply given by the final two terms of Eq. (48). The reason is that in 
this case the solution of the generalized GPE has a different shape from the solution of the ordinary GPE. Thus the 
condensate interaction terms Afo(00|y |00) are numerically different in the two equations, even though they have the 
same functional form. In the homogeneous limit this issue does not arise and we can write an explicit expression for 



AEx{p) [see Sec. VI B|. We note also that the energy shift AEx{p) arises naturally as a consequence of the constraint 
on N which means that an excitation can only be created by removing particles from the condensate. However, the 
same shift also appears in a broken symmetry (and hence grand canonical) approach because the chemical potential 
is taken to be Ag rather than A for the higher order calculation. 

The energy shifts described by AE^ij)) and AE\(jj) are quadratic forms in the quasiparticle transformation coeffi- 
cients Upi and Vpi for the particular mode p under consideration [this follows from the fact that Apjiijp) and Anjiij)) 
are quadratic forms, c.f. Eq. (|3l|)]. They can therefore be written straightforwardly as a modification to the matrices 
C and which appear in the BdG equations of Eq. (H) [c.f. Eqs. (H), (|l|) and (||)]. The new matrices are defined 

by 

^u" = Htf ~ XgS,, + 2No{0i\V'm + E 2{kt\V'\jm)pmk, (65) 
= No{ij\V'm + J2 {^j\V'\km)Kk„,. (66) 

fcm/O 



Writing the energy shifts in the form of Eqs. (|6g) and (|6g) allows the perturbative calculation to be made self- 
consistent if an exact diagonalization of the new BdG equations is performed. When we refer to ordinary perturbation 
theory we will therefore mean the simple evaluation of the expressions for Ai?4(p) and AE\{p) of Eqs. ( pi] ) and (|6^ ) 
using quasiparticle energies and transformation coefficients calculated from the quadratic BdG equations [i.e. with £ 
and Ai given by Eq. (^8|)]. Self-consistent perturbation theory, on the other hand, refers to the exact diagonalization 
of the generalized BdG equations which are obtained from the matrices C and A4 of Eqs. ( |65| ) and (^6|) (and also 
including the corrections from ^3 which will be discussed in the next section). 

Comparison of Eqs. (p^), ( |65| ) and ( |66| ) shows that there are three types of correction introduced by Ai?4(p) 
and AE\{p). First, the condensate energy A which appears in C is upgraded to the value appropriate to the 
generalized GPE. Second, direct and exchange collisions between noncondensate atoms are included via the term 
^j.j^_^g 2(fcz|y |jm)pmfe which appears in £y . Finally, the anomalous average k^- is introduced into the off-diagonal 
terms Mij. We will show in Sec. ^ that the interpretation of this correction is that it upgrades the condensate- 
condensate interactions which appear in the leading order contribution to Mij so that these are described by a 
many-body T-matrix. 

We note that the changes to the coefficients Cij and A4ij which are introduced by Ai?4(p) are exactly the same 
as would have been obtained by using the factorization approximation of Eq. ( |5^ ) on the operators appearing in 
H4. Thus this approximation on the product of four operators is equivalent to self-consistent first order perturbation 
theory and is therefore justified to the order of this calculation. This result is not too surprising given that Wick's 
theorem is concerned with operator averages and hence is closely related to first order perturbation theory. 

Although the expressions for Ai?4(p) and AE\{p) have been derived for pure quasiparticle number states, they can 



be straightforwardly reinterpreted as thermal averages along the lines of the discussion in Sec. HI D. The energy shifts 
therefore become functions of temperature via their dependence on the quasiparticle populations and transformation 
coefficients. 
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D. Second Order Perturbation Theory 



The expression for the energy shift of a quasiparticle number state \s) from second order perturbation theory is 

i;Pcrt(s, 2) = 2^ ^ _^ , (67) 

r^s ^ ^ 

where Eg and Er are energies of the system calculated from the quadratic Hamiltonian via Eq. (|33|). As mentioned 
ear her, we can neglect the contribution from H4 at this order of perturbation theory and consider only the modified 
cubic Hamiltonian defined by iJg — H3 + AHi . For this reason we will denote the second order shift in the energy as 
E3. 



It is convenient to rewrite the Hamiltonian in the quasiparticle basis. Using Eqs. (pi]), (21) and ( p^ we obtain 
the result 



where the coefficients Aijk, Bijk and C'i are given by 



+ h.c. , (68) 



(69) 



7nnqy^0 

+ {mn\V''\qO) [U*,U,,Vkrn + U*„,V,nUkg + V,lV,nVkm] , (70) 

Ci + ^(^(jgi + ^(jiq)"-(; = 0. 
q=iO 

These coefficients can be rewritten straightforwardly in the position representation [where they have the form of 



integrals over the quasiparticle functions u{r) and v(r), c.f. Eqs. (|8l|)-(83)] using Eq. ([36[). We have given the linear 
coefficients Ci in the particular combination of Eq. (Flj) because they appear in the energy in this form [see Eq. ([72|)]. 
Linear terms in the Hamiltonian arise both from AHi and from the use ofquasiparticle commutation relations when 
H3 is written in normal order. The fact that the right hand side of Eq. (|7l| ) is zero demonstrates that AHi cancels 
naturally with part of the contribution to the energy from H3. In fact, AHi removes that part of which is obtained 
from the factorization approximation as we will discuss near the end of this section. 
The modified cubic Hamiltonian gives a contribution to the energy which is 



E: 



p 12 



3 



1 I Ji I 

7 [ntUjUk - {rii + l)[n.j + l){nk + 1)] 

6 .i^o + + 

^ 1 ^ \B.,k^B,k,? ^(^^ ^ ^^^^^^ _ ^^^^^ ^ ^^^^^ ^ 

- E ^IC*'' + Y^B,,^ + Bnj)n,\\ (72) 

where we have again left out terms which are negligible in the thermodynamic limit (the full expression is given in 
Rcf [HI). The coefficient A^-^ is defined as a sum over the permutations of the three indices in Aijk, i.e. Af-^. = 
Aijk + Ai]^j + Ajik + Ajki + Akij + Akji- We note that the expression for E^ contains quasiparticle population terms 
in exactly the factors we would expect for the Bose enhancement of scattering into and out of the states involved. 

Equation ( [72[ ) gives the contribution to the total energy of the system from second order perturbation theory. The 
contribution to the energy of a quasiparticle corresponds to the change in E^ when a quasiparticle is created in some 
particular mode p. We are therefore interested in the change in £'3 as np + 1, which we denote by AEs{p) 
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AEs{p) = E3{ni,n2 ...Up + l...)- E3{ni,n2 . 
Evaluating AE^lp) using Eq. (p2|) gives 



(73) 



AE,ip) = - ^ 



2(ep- 






+ Bpji 1 


2(ep- 


- e, - cj) 


1 ^ijp 


+ -Sipj 1 




e» + ej 



[1 + rtj + Uj] 



[1 + + flj] 



(74) 



Although this expression has been derived for pure states, we will now reinterpret it as a thermal average along 
the lines of the discussion in Sec. HID. As before, this has the effect that the energy shift becomes dependent on 
temperature. 

We will show in Sec. ^ that the physical interpretation of Ai?3(p) is that it introduces a T-matrix into the description 
of condensate-noncondensate collisions in the diagonal terms of the BdG equations. We note here, however, that the 
first term in Eq. ( |7^ ) corresponds to the simultaneous annihilation or creation of three quasiparticles, while the 
second describes Beliaev processes in which a single quasiparticle spontaneously decomposes into two others 
These processes can occur at zero temperature and are therefore dominant in the low-temperature regime. The final 
term corresponds to Landau processes in which two quasiparticles collide and coalesce to form a single quasiparticle 
(these are essentially the reverse of Beliaev processes). Landau processes can not occur at zero temperature (because 
there are no excited quasiparticles) but they dominate at high temperature. If an energy matching condition is 
satisfied, then the Beliaev and Landau processes can occur in a real rather than a virtual sense and lead to a finite 
lifetime for the quasiparticles. In the homogeneous limit these lifetimes can be calculated analytically (Sec. VIB 2| ) 
and give the same result as the application of Fermi's Golden Rule to the Hamiltonian H^. In a trap the situation is 
more complicated because there may be no exact energy matches, and only a few states which are close enough to be 
strongly coupled. Nonetheless, A£'3(p) can be used for a trapped gas to calculate the time evolution of a quasiparticle 
and determine its lifetime if this is meaningful. Such calculations are currently underway and will be described in a 
forthcoming publication 139]. 

The expression for AE^^p) of Eq. (|7j) is a quadratic form in the quasiparticle transformation coefficients Upi and Vpi 
for the particular mode p under consideration. This can be seen from the fact that each of the indices in the coefficients 
Ap^j and Bpij corresponds to a linear dependence on a quasiparticle transformation matrix U or V. Although Eq. ( [T^ ) 
is all that is required for ordinary perturbation theory, if we wish to calculate Ai?3(p) self-consistently we need to 
rewrite it in a form which explicitly separates out the quadratic dependence on Upi and Vpi. The result is 



AE^ip) 



AC{ep) AM{ep) 
-AM*i-ep) -AC*i-ep 



(75) 



where the matrices AC{ep) and AA4{ep) depend on the energy of the quasiparticle mode ep (and on temperature) 
and are defined by 



km^O 



AM,,iep)= J2 

km^O 



(1 + rifc + Hm) 



jkm 



— {^k + ^m) + + 



(76) 



(1 + rifc + Um) 



AikmBjkrri ^*fcm^jfcn 



Cp — (ffc + Cm) Cp + Cfc + Cm 
, ("-f" ~ "-fc) ^ikm^jmk 



^p ~ 



where 



qr^O 



(77) 



(78) 
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B^km = {<l^\y"\^0)VkqVmr + {rO\V'\qi)[UkgVrnr + VkrUraq] 



qr^O 



qr^O 



(79) 
(80) 



Aikm = 2y/l%Ua J d^r C*(i")|Co(r)"fc(i")'«m(r) + Co(r) Ufe(r)w„(r) + Vk{r)um{r) 
Bikm = 2^/NoUQ J d^r Q{r)l^(Q{r)vk{r)vm{r) + Co(i") Mfc(r)w™(r) + Vk{r)ujn{r) 
Qfcm = 2y^Ua I dhQ{r)\Ca{r)ul{r)v„,{r)+C*{r)\ul{r)u„^{r)+vl{r)v„^{r) 



These coefficients can be written in tlie position representation using the contact potential (in which form they may 
be more convenient for calculations) as 



(81) 
(82) 
(83) 



The quantities A£y(ep) and AMij{ep) satisfy A£*^{ep) ~ A£y(ep) and AA4ji{ep) — AMij{~ep). 

Since AE3{p) simply modifies the matrices C and A4 which appear in the BdG equations, its calculation can be 
made self-consistent by including A£(ep) and AA4{ep) in the BdG equations and solving them exactly. If we also 
include the effects of AE4 [p] and AE\ (p) self-consistently, then the generalized BdG equations are 



^ep) 



Miep) 



where C{ep) and M{ep) are the matrices with elements 

C^j{ep) = C^J + AAj(ep), -^ij(ep) = + AA^y (ep), 



(84) 



(85) 



and Cij and A4ij 



are defined in Eqs. ( pSD and {p6\j, while ACij{ep) and AMij{ep) are defined above. Eq. ( p4D defines 
the final form of the BdG equations to the order of the calculation presented in this paper. 

We note that the inclusion of AE^lp) modifies the structure of the BdG equations, as can be seen from a comparison 
of Eqs. ( p2|) and j^^j). There are two important changes. The first is that the matrix which is to be diagonalized now 
depends on the energy of the quasiparticle mode under consideration. This means that a single matrix diagonalization 
no longer yields the whole quasiparticle spectrum. The second change is that (for 7^ 0) the diagonal elements are no 
longer proportional to (the complex conjugates of) each other, and similarly for the off-diagonal elements. However, 
this proportionality is a consequence of any quadratic Hamiltonian, i.e. any Hamiltonian of the form of Eq. (fTT 



with arbitrary coefficients Cij and A4ij. Thus the full effect of Ai?3(p) can not be reproduced by any quadratic 
Hamiltonian. This is the reason why a variational approach to the problem of the dilute Bose gas only reproduces 
the HFB theory and does not include the effect of A£'3(p) (see Sec. |IVED [0. 



A further feature of AE^^p) is that it is intrinsically non-local. Thus even if a contact potential is used to describe 
particle interactions, the position representation of Eq. (|7^) still involves integrals over two spatial coordinates. In 
contrast, the energies obtained from the quadratic theory or from Ai?4(p) depend only on a single spatial integral if 
the contact potential is us ed. 

As we discussed in Sec. EVB, the conventional treatment of the cubic Hamiltonian is based on the factorization 
approximation. A comparison of Eqs. (0), (|5^ ) and ( ^3|) shows that this leads to a modified cubic Hamiltonian 
— H3 + AHi which is identically zero, so E3 and Ai?3(p) are completely neglected. This means that the 
factorization approximation only takes into account that part of which cancels with AHi. This term appears 
as a result of the change in shape of the condensate when we use the generalized rather than the ordinary GPE. 
Thus the factorization approximation takes condensate shape effects into account but neglects all the Beliaev and 
Landau processes which can occur in the noncondensate. These give a contribution to the energy which is of the same 
order of magnitude as AEi{p), however, so the factorization approximation results in an inconsistent treatment of the 
non-quadratic Hamiltonian. Indeed we will show in Sec. that AE3{p) is required to remove infra-red divergences 
in the theory and leads to a gapless spectrum. The factorization approximation on the product of three operators is 
therefore invalid and is the cause of many of the difficulties encountered in extending the Bogoliubov theory of BEC 
to higher order. 
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E. HFB and HFB-Popov 



The generaUzed HFB theory includes the effect of the non-quadratic Hamiltonian via the factorization approxima- 
tions of Eqs. (|5^) and (^4|). As we have seen, these approximations correctly include the contribution from AE4 but 
neglect the terms in AE3. The generalized HFB theory therefore involves solving the BdG equations in the form of 



Eq. (22) but with the matrices C and A4 taken from Eqs. ( pq ) and (66). The condensate wave function and eigenvalue 
are determined from the generalized GPE of Eq. (^ ) . 

The HFB theory can also be obtained using a variational approach . In this case the Hamiltonian is assumed to 
have the quadratic form of Eq. (|l^), but the coefficients Cij and A4ij are treated as variational parameters which are 
used to minimize the free energy of the system. This procedure leads to the same coefficients as are obtained by the 
factorization approximations and hence this approach reproduces the HFB theory. The variational calculation does 



not include the effect of AEsijj) because, as we mentioned in Sec. iVD, this leads to an expression for the energy of 
a quasiparticle which does not correspond to any quadratic Hamiltonian. 

The full HFB theory is not used in practice, however, because it does not predict a gapless spectrum and it suffers 
from infra-red divergences (see Sec. |Vl| ). Instead, the Popov approximation is used in which the contribution from the 
anomalous average is neglected in both the generalized GPE and the BdG equations. In the position representation 
(using the contact potential), the BdG equations therefore have the form of Eq. js^ ) with C{r) and A^(r) given bjj^ 

£(r) = H'P - Ag + 2[/o[7Vo|Co(r)P + Pe.(r)], M{r) = NoUoC^r). (86) 



The HFB-Popov theory is gapless and it is also free of infra-red divergences as we will show in Sec. VI. It is therefore 
preferable to the full HFB approach and for this reason it is the basis of recent numerical calculations at finite 
temperature [ ^Opl| ]. 

An interesting feature of the Popov approximation is that the shifts in the quasiparticle energies depend predom- 
inately on the spatial variation of the noncondensate density in the region of the condensate. This can easily be 
seen if we write the expression for AE\{p) of Eq. (^) in the position representation using the contact potential 
approximation. If we ignore the fact that the condensate shape changes in going from the ordinary to the generalized 
GPE, then AE\{p) is given by 

AEx{p)^-Uo jd'v2p,,X^)\Uv)? jd'v'Appiv'). (87) 



Comparison with the expression for Ai?4(p) of Eq. ( |62| ) (with K(r) 0) shows that if Pex{f^) is independent of r then 
the combination AE4^{p) AE\{p) is zero. A consequence of this is th at in the homogeneous limit there is no change 
in the excitation spectrum within the Popov approximation [c.f. Eq. (n2q)]. 



V. ULTRA-VIOLET DIVERGENCES AND THE T-MATRIX 

In the previous sections we have used the bare interatomic potential V(r) to describe particle interactions. It is 
well-known, however, that at low temperatures the scattering of neutral atoms in a three dimensional gas can be 
characterized by the s-wave scattering length a. This parameter is usually introduced by replacing V{r) with the 
contact potential of Eq. (|I|) , although this leads to the appearance of ultra-violet divergences in the theory. This is 
not surprising when we consider that a contact potential can scatter high-energy atoms as effectively as low-energy 
ones. Of course this is physically unrealistic, and in reality the momentum transfer between atoms will vanish at 
large momenta {k > 1/a). The contact potential is a low-energy approximation and care must be taken to ensure 
that high-energy states are dealt with correctly. 

The purpose of this section is to enable the contact potential to be used to describe particle interactions while 
introducing the required ultra-violet renormalization in a rigorous manner. This is achieved using the fact that the 
contact potential is really an approximation to the low-energy limit of the two-body T-matrix (T2b) which describes 
particle scattering in a vacuum. We will therefore rewrite interaction matrix elements in terms of this T-matrix 
and show that the difference between and the interatomic potential ^(r) provides the necessary ultra-violet 
renormalization. 



®The equations are usually given without the coefficients Cj, i.e. as in Eq. (^l[). Excitations orthogonal to the condensate can 
be obtained in this case using the prescription of Eq. (^). 
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We note, however, that the procedure we will describe is of general validity and can be used for situations where 
the contact potential approximation to the T-matrix does not apply. In particular, the expressions we will obtain can 
be used to describe a charged Bose gas or neutral atoms in two dimensions and the only change is that a different 
substitution for the T-matrix is required. For this reason we have written all equations in terms of T2h, and only 
replaced this with Uq S{r) when we have specialized to the case of neutral atoms in three dimensions. 

As well as the removal of ultra-violet divergences, there are a number of additional reasons why it is important 
to rewrite matrix elements in terms of T2h- The first is the fact that the details of the interatomic potential are 
often not very well-known, while the low-energy T-matrix is fairly universal and can be characterized by a single 
parameter, namely the s-wave scattering length a. The second is that the contact potential is very convenient for 
numerical calculations since it leads to a considerable simplification of the equations. A third reason is that, for 
singular potentials (such as a hard sphere), the matrix elements (ijlV^lkm) are actually very large or poorly defined 
for low-energy states. This means that the Hamiltonian written in terms of V{r) is not convenient for numerical 
calculations. In contrast, the T-matrix elements describing scattering off such potentials are finite and well-defined. 

Perhaps the most important reason for the introduction of the T-matrix, however, is the fact that our perturbative 
analysis can not be expected to be valid if interactions are characterized by the actual interatomic potential V{r). The 
reason is that the physical interpretation of higher order terms [specifically Kij and AE^p)] is that they introduce 
T-matrix corrections into the description of particle collisions, as we will show in Sec. |V C| . However, a perturbative 
treatment of two-body collisions (i.e. the Born approximation) is known to fail at low-energy. This is apparent from 
the fact that this approach gives ~ ^(r), whereas the low-energy T-matrix is actually independent of the details 
of V{r) (as in the contact potential for example). It is therefore necessary to take into account those terms which 
upgrade V{r) — > to all orders and regroup terms so that the Hamiltonian is written in terms of T2h- We will show 
that perturbation theory then corresponds to a Born expansion of the difference between scattering in the presence 
of a condensate (described by the many-body T-matrix T„ib) and scattering in a vacuum. These effects are small for 
a dilute gas away from the critical point and therefore a perturbative treatment should be valid. 

The structure of this section is as follows: in Sees. V A and V B, we give a brief discussion of the two-body and many- 
body T -mat rices. In Sec. V C we discuss the physical interpretation of Ktj and AE3{p) in terms of these quantities. 



In Sec. VD we show that the Hamiltonian can be rewritten in terms of T2h, and that this leads automatically to an 
ultra-violet renormalization of the theory. Some of the mathematical details of the argument are given in Appendix 



A. The Two-Body T-Matrix 

The two-body T-matrix is defined as a function of the complex parameter z by the Lippmann-Schwinger equation 

T2b(^) - ^ + E ^N)'" , si ^ sp, ^''(TO|T2b(^), (88) 

pq ^ ^'^P + ^1 > 

where and e^^ are particle energies (eigenvalues of H^^). The kets \pqY^ are the corresponding two particle 
eigenstates and describe the intermediate states in the collision of two atoms. The summation over these states 
includes the noninteracting ground state.P^ 

In the homogeneous limit, the scattering between states with relative momentum ?ik and ?ik' is independent of the 
centre-of-mass momentum KK and is described by the matrix element r2b(k',k, z) = (k' K|T2b(z)|kK). Physical 
scattering events correspond to on-shell matrix elements for which k = k' and z — lim^^o + i-V (here k = |k|, 
e^^ = h^k^ /2ifi, rh = m/2 is the reduced mass, and we have redefined z to include the centre-of-mass energy). For 
neutral atoms in three dimensions, these matrix elements can be calculated in the low-energy limit fca ^ 1, with 
the result liiii|i(;|=|k'| '^2b(k', k, e^^ -I- irj) = Uq + O(fca). This implies that the spatial representation of the low-energy, 

on-shell T-matrix is given by the contact potential of Eq. (|1|) . Measurements of the scattering length thus correspond 
to measurements of this T-matrix. 

We note that the contact potential is only a valid approximation at low-energy and the full T-matrix vanishes for 
k > 1/a. The simplest way to include this in the theory is to use the contact potential together with a momentum 



^"The labels of the intermediate states in Eq. ( |88| ) correspond to single-particle wave functions (^^''{r). For a trapped gas these 
are different from the basis functions (i{r) we have used in the theory so far because they are not orthogonal to the condensate. 
This does not create any difficulties, however, since Eq. (^) simply serves to define an operator whose matrix elements become 
an input to the theory. 



19 



space cut-ofF around k = 1/a. In fact our results do not dep end on the position of this cut-off so that it may ultimately 
be taken to infinity. This result (which is discussed in Sec. [VD| ) demonstrates that there is no ultra-violet divergence 
in the theory if the contact potential is used as an approximation to T2h rather than ^(r). 

Although the relationship between the two-body T-matrix and the contact potential has only been discussed for 
on-shell elements in the homogeneous limit, we will assume that the result can also be used to calculate off-shell 
matrix elements, and in addition that it can be applied to a trapped gas. In the latter case, this assumption relies 
on the validity of a local density approximation for two-body scattering in a trap and is discussed in more detail in 
Rcf. H. 



B. The Many-Body T-Matrix 

In a Bose condensed system, particle collisions occur in the presence of the condensate and are therefore described 
by a many-body T-matrix. In this section we will define this T-matrix, discuss the physical effects it describes and 
relate it to the two-body T-matrix. 

The many-body T-matrix T^\j{z) appropriate to this paper is defined by the Lippmann-Schwinger equation 

where z is a complex parameter^] and €p and Hp are quasiparticle energies and populations. We note, however, that 
the intermediate states \pq)'^^ correspond to particle wave functions (see below). In the homogeneous limit this does 
not pose any difficulties because a single momentum can label both particle and quasiparticle states. For a trapped 
gas, however, the definition of Eq. (|8^) only applies in the high-energy limit where quasiparticle and particle wave 
functions are identical. 

Comparison of Eqs. ( |8^ ) and (|8^) shows that the many-body T-matrix takes into account two effects of the medium 
in which collisions occur. First, the intermediate states may be occupied so there is Bose enhancement of scattering 
into and out of these levels described by the factor 1 + Up + nq = (rip + l)(np -|- 1) — UpUq. Second, the existence of 
a condensate means that the energies and populations of the intermediate states correspond to quasiparticles rather 
than particles. 

The full quasiparticle character of the intermediate states is not taken into account in the T-matrix of Eq. 
however, because the transformation coefficients Uij and Vij do not explicitly appear. Thus the wave functions of 
the intermediate states are treated in the perturbative (particle) approximation Uij 5ij, Vij ^ 0. A more general 
T-matrix which includes the full quasiparticle wave functions has been discussed by Bijlsma and Stoof |2|]. It is the 
simpler expression of Eq. ( |89| ) which appears naturally in the theory, however, at least to the order of this calculation 
[c.f. Eq. (11)]. ^ ^ 

It is convenient to rewrite T^b in terms of T2h because both these quantities are well-defined for singular potentials. 
Using their respective definitions, it is straightforward to show that they are related by 

rr ( \ rr. i-^.sr^ J'2b ( ^) | w) ( 1 + Tip + Tig) {pq\T^,h{z) 

~J2T2bi-z)\pqr - fsl^ sp, ^^(Wlrmb(^). (90) 
pq ^ ~ 'y^P + ^9 

This result greatly simplifies the calculation of T,„b in terms of the s-wave scattering length. In addition, although a 
Born expansion of T2h and Tmb in terms of V{r) will fail at low energies, an expansion of Tmb in terms of T2b based 
on Eq. (po[) should be a good approximation for a dilute gas away from the critical point. 



Although the many-body T-matrix is defined for a general complex z, this parameter is physically interpreted as the energy 
of a collision. We note therefore that the zero of energy in Eq. (p^) is not the same as in Eq. (fes ) . The many-body T-matrix 
is defined in terms of quasiparticle energies so z is measured relative to the condensate in Eq. (p9). The two-body T-matrix is 
written in terms of particle energies so in that case z is measured relative to the energy of a stationary particle at the bottom 
of the trapping potential. 
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C. Interpretation of Kij and AE3 

In this section we will discuss the physical interpretation of and AE^ij)) and show that they introduce the 
many-body T-matrix defined above into the description of particle collisions. Our discussion is based predominately 
on the homogeneous limit, although the results also apply to a trapped gas at high energy. 



1. Interpretation of itij 

The anomalous average is defined in terms of the quasiparticle transformation coefficients by Eq. (|3^), and as 
far as the numerical implementation of the theory is concerned it is simply calculated directly from this definition. 
For the purposes of interpretation, however, it is useful to express it in terms of quantities which have a more direct 
physical significance such as interaction matrix elements and quasiparticle energies and populations. This can be 
achieved in the high-energy limit 00) using Eq. ( p5| ) which gives 

{l + n, + nj)M^j 

^'^^ -(e. + e,) ' ^'^^ 

In the quadratic theory, is given by Eq. (|8|), so comparison with Eqs. (^ ) and ( |66|) shows that Kij introduces 
the second order Born approximation to the many-body T-matrix into condensate-condensate interactions. 

In fact, self-consistent perturbation theory leads to the introduction of this T-matrix to all orders in its Lippmann- 
Schwinger definition. This result can be demonstrated in the homogeneous limit where the BdG equations which 
determine Hij can be solved analytically. In this case it is easy to show that Eq. is exact for all provided 
only that Mij does not depend explicitly on ep (which is the case if we neglect /S.E^{p) or approximate it by setting 
Cp 0). In particular, the result applies if Ai?4 is calculated self-consistently, in which case Mi^ is given by Eq. (|6^). 
This leads directly to the results 

M.,j - A^o(«j|r,„b(2 = 0)|00), (92) 



= ^^^^i^^fi^iVo(*j|r,„b(^ = 0)|00), (93) 

with Tmb(2:) defined as in Eq. (|89|). 

The interpretation of the anomalous average is therefore that it introduces the many-body T-matrix into the 
description of condensate-condensate collisions in both the generalized GPE and the off-diagonal elements of the 
BdG equations. Self-consistent perturbation theory leads to the introduction of this T-matrix to all orders in the 
Lippmann-Schwinger equation, while ordinary perturbation theory introduces the second order contribution in a Born 
expansion. 



2. Interpretation of AEs 



Although the general expression for AE^ is rather complicated, the physical interpretation of this term can be 



seen by considering its form when the coefficients Aijk, Bijk and Cijk of Eqs. ([ 

(high-energy or perturbative) approximation [/^ Sij, Vij 0. In this case Eqs. (76|) and ( |77| ) become 



AC^j{ep) 



E 

km^O 



2No{Oi\V''\km){l + rik + n^){km\V^\jO) 



{f-k + Em) 



|) are treated in the particle 
,nd (|7^) become 

mo{im\V'\kQ){n^ - nk){kQ\V'\jm) 



Cp -\- Cn 



ek 



and 



AM.jiep)^ J2 

km^O 



ANoiimlVmin^ - nk){jk\V'\mO) 



(94) 



(95) 



shows that the first term in ACij(ep) is the second order contribution to Tmb for the 



Comparison with Eq. 

expression 2iVo(0i|y |jO), which describes condensate-excited state collisions in Cij [c.f. Eq. ([18|)]. The T-matrix is 
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introduced with the on-shell value of z since in the limit e^, — s- cx) the quasiparticle energy is given by Cpp [c.f. 

Eq. (§)]. 

The second term in the equation for A£ij(ep) is similar in structure to AA^ij(ep). Each of these expressions 
contains the effect of four distinct scattering processes because they are written in terms of symmetrized matrix 
elements. The simplest contribution to l^Mij{ep) is illustrated in Fig |], and is a modification to the processes of 
pair production out of the condensate described by the leading order expression for M.ij of Eq. ([l8|). Figure [l] shows 
that AA^ij (ep) describes interactions between condensate atoms which occur via the exchange of a particle from the 
noncondensate. These interactions can also be thought of in terms of polarization of the medium in which collisions 
occur; the first condensate atom causes a density fluctuation in the noncondensate which then affects the motion of 
a second condensate atom. 

Thus, just as introduces the many-body T-matrix in the off-diagonal terms of the BdG equations, Ai?3(p) 
introduces these corrections into the diagonal terms and in addition introduces more complicated polarization effects. 
An important difference, however, is that whereas self-consistent perturbation theory on introduces Tmb to all 
orders, Ai?3 only introduces it to second order in a Born expansion. This result is not altered by making the 
perturbation theory self-consistent, because the result of Eq. (^4|) was derived by considering the high-energy limit 
Uij Sij, Vij 0. This limit is independent of the form of the interaction potential and is therefore unchanged 
as the perturbation theory is made self-consistent. The terms which introduce the third order contributions to T^h 
in condensate-excited state collisions are in fact beyond the order of this calculation and come from the use of third 
order perturbation theory. The fact that AE^ only introduces T^b to second order has implications for the gapless 
nature of the theory if self-consistent perturbation theory is used and is discussed further in Sec. VI. 

The anomalous average and AE^ are the terms in the GPE and BdG equations which contain ultra-violet divergences 
if the contact potential is used as an approximation for V(r). The fact that both these quantities introduce T-matrix 
corrections suggests that there is an intimate connection between ultra-violet divergences and the T-matrix and 
provides a simple explanation of why such divergences occur. As we discussed in Sec. V A, the contact potential is 
really an approximation to T2h rather than to V(r). Some of the corrections introduced by Kij and AE3 exist even 
in the extremely dilute limit and correspond to the introduction of T2h- It is therefore incorrect to use the contact 
potential together with the full expressions for Kij and AE^ as this amounts to double counting two-body effects. We 
will show in the next section that the correct use of the contact potential as an approximation for T2h avoids this 
double counting and leads directly to a renormalization of both Ky and AE^ . 



D. Ultra- Violet Renormalization 



In this section we show how all interaction matrix elements can be rewritten in terms of the two-body T-matrix 
and give explicit expressions for the ultra-violet renormalization of and AE3 (p) . The mathematical details of the 
argument are given in Appendix |^. 

The introduction of the T-matrix occurs in two stages. In the first, we return to the original Hamiltonian of Eq. (j^) 
and divide the single-particle basis states Ci(r) into two subspaces, a low-energy region (L) and a high-energy region 
(H) . In the low-energy subspace interactions between particles may be significant and the quasiparticle techniques we 
have discussed in the previous sections are needed. In the high-energy subspace, however, the effect of interactions is 
small and the characteristic time scale of evolution is fast compared to the low-energy states. It is therefore possible 
to integrate the equations of motion for high-energy operators and express them in terms of operators which act only 
on low-lying states. This results in the replacement of the bare interaction potential between two low-energy particles 
with an effective interaction which has the form of a restricted two-body T-matrix, the restriction being that the 
intermediate states must be of high energy. 

T he iri itial Hamiltonian of Eq. is therefore replaced with an effective low-energy Hamiltonian which is [c.f. 
Eq. ( [ATI )] 



L ^ L 

HeS = ^ H-^alaj + 2 X! I^h| ^"1)0! a] afca„, (96) 

ij ijkm 

where means a sum over low indices. The restricted tw o- bod y T-matrix Tjj which describes low-energy interactions 



is defined by the Lippmann-Schwinger equation [c.f. Eq. (|A15| ) 



H 



Th^V + Y, Vlpqr" , ^ sp, ''{pq\TH, (97) 



22 



which is the usual expression for the two-body T-matrix [c.f. Eq. (88)] except for the fact that the intermediate 



states are restricted to the high-energy subspace. All our previous results can therefore be taken over by making the 
substitutions {ij\V\km) (ijlTulkm) and J2 ^ (matrix elements involving Th are to be interpreted as being 
symmetrized). Since the Hamiltonian now has an explicit high-energy cut-off, the problem of ultra-violet divergences 
is replaced by the problem that physical quantities may depend on the cut-off frequency. 

The second stage of the procedure is to replace the restricted T-matrix Th with the full two-bo dy T -matrix of 



Eq. (Bq). This is achieved using the fact that these quantities are related by the expression [c.f. Eq. (Al£)] 



T2b{z) =Tu + Y, Talpgy , sp ^ sp. '"{pqlT^biz). (98) 

pq \ f H / 

The importance of this result is that Th should be approximately equal to T2h for systems where the density and 
temperature are low enough that the boundary between high and low-energy states can be taken fairly low. In this 
case the second term in Eq. (^8|) should be small compared with the first, and a Born expansion of Th in terms of T2h 
should be valid. This is in contrast to the fact that a Born expansion of either Th or T2h in terms of V^(r) is expected 
to fail at low energy and is essentially the reason for the introduction of ThJ^ We show in Appendix]^, that at T = 
the Born expansion of Eq. ( p8| ) is valid in three dimensions if the dilute gas criterion na^ ^ 1 is satisfied. At finite 
temperature the requirement becomes a/XdB ^ 1 where XdB is the de Broglie wave length. This condition is satisfied 
for a dilute gas even at temperatures above the critical point. 

We therefore introduce the full two-body T-matrix by writing Eq. (^3) in a second order Born approximation as 



Th ~ T2i,{z) - VT2b(z)|M)'^ , , sp, {pq\T2biz). (99) 



\sp 

pq 



The contribution from the first term amounts to replacing Th with T2h{z) in all interaction matrix elements. This is 
all that is required (to the order of this calculation) in and H4, but in Hq and H2 the second term must also be 
included. This has the effect of renormalizing both and Ai?3(p) so that they are finite even if a contact potential 
is used for T2b- The rcnormalization therefore simply amounts to the addition and subtraction of the second term 
in Eq. (p9|). The addition of this contribution gives a better approximation to T2b in leading order matrix elements, 
while the subtraction renormalizes higher order terms and ensures that two-body effects are not double counted. Since 
the theory is now well-behaved in the high-energy limit, all summations are convergent and they may therefore be 
extended to infinity. This means that physical quantities do not depend on the exact position chosen for the cut-off 
between high and low-energy states. 

Although this procedure is valid mathematically for any (small) value of z, the physically appropriate values are 
those which arc introduced naturally by and AE3{p) in the two-body perturbative limit. In this case we have 
Ep — [c.f. Eq s. ([l3| ) and (^], so Eq. ( |9l| ) shows that in the matrix elements {ij\T2h{z)\00) we should 

set z = 2eQ^, while Eq. ( |94| ) shows that for the matrix elements {Oi\T2h{z)\jO) we should set z — e^P + e^^ . In the 

highest order matrix elements (those introduced by H3 and H4), the values of z should be taken to be the two-body 
limit of the energy of the corresponding collision process. This means that in terms involving k^- we have z = 2eg'', 
because k^- always introduces the T-matrix into collisions of the form (ij|T2b(2;)|00). Terms involving p,y are more 
complicated and involve a range of different energies. The appropriate values of z can be found by decomposing 
into quasiparticles via Eq. (|29|), and considering the two-body limit of the various different quasiparticle collisions 
which then appear. In three dimensions, however, the T-matrix is independent of z at low energy and is simply given 
by the contact potential. 

We have now succeeded in rewriting all interaction matrix elements in terms of the two-body T-matrix. For neutral 
atoms in three dimensions this can be replaced by the contact potential in the honiogeneous limit. For a trapped gas 
the situation is more complicated, however, because the T-matrix defined in Eq. (|8g) is not the same in a trap as it 
is in the homogeneous limit. Nonetheless it is this T-matrix which appears in the theory so in principle we should 
calculate all matrix elements of T2biz) explicitly from their definition. This is hardly practical, however, and we will 



^^Difficulties will arise, however, if T2\^{z) depends strongly on z at low energy, or (equivalently) if Th is sensitive to the 
position of the cut-off. In this case the Born approximation [which in leading order is simply Th ~ T2h{z)] will fail. This is not 
a problem in three dimensions where T2]:,{z) tends to a constant (the contact potential) at low energy, but it will be important 
in two dimensions where T2y,{z) ~ — l/ln(z) as 2 ^o|. In this case we will have to work explicitly with Th. 
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assume instead that the T-matrix in a trap can also be replaced by a contact potential. This assumption relies on 
the validity of the local density approximation and should be valid if the dominant contribution to particle scattering 
involves intermediate states of high energy which are insensitive to the details of the trapping potential. This issue is 
discussed further in Ref. Q. 



1, The GPE and BdG equations 
Following the procedure discussed above, the generalized GPE of Eq. ( ^ ) becomes 

Hll + 7Vo(fcO|T2b(2eg^)|00) + ^ [2(fcz|T2b(z)|jO)p,, + (fcO|T2b(2eg^)|ij>j',] = Ac 4o, (100) 

where the values of z in the term involving pji are calculated as discussed above. The renormalized anomalous average 
nfj is defined by 



where is the two-body, perturbative limit of [c.f. Eq. (|9l|)] 

J^2b(2eg^)| 



ij ~ ^sp isp . sp\ ■ V-*^"^; 



In the position representation using the contact potential, Eq. (101) becomes 

K^(r) =K(r)-«;"P(r), (103) 

with 

The renormalized form of the coefficients Cij{ep) and Aiij{ep) of Eq. ( |85| ) is 

Aj(ep) = H^^ - AG<5y + 27Vo(0^|^2b(6^^ + e^f )|jO) + ^ 2(fci|T2b(z)| jm)p,„fe + A£f^{ep), (105) 

X.,(ep) = 7Vo(zj|T2b(2eg^)|00) + J] (zj|r2b(2e;;^)|fcm)<™ + ^M^A^p), (106) 

where A£,^ (ep) is calculated from the renormali zed f orm of AE^{p) (see below) and the condensate eigenvalue Xq 
is calculated from the renormalized GPE of Eq. (IOC). The values of z in the term involving pmk are calculated as 
discussed above (Sec. VP ). The renormalization of A£'3(p) comes from the use of the second order expression of 
Eq. ( p^ ) in the matrix element 2A'^o(Oi|T2b(eo^ + ep^)b0) which appears in £y (tp). The new form of AE3{p) is 

AE^ip) = Ai?3(r2b,p) - Ai?f (p), (107) 

where 

^E/{p) = 2No^^ — ,p Apjiip). (108) 



Here Ai?3(T2b,p) stands for the expression of Eq. ( |74[ ) with the coefficients Aij^ and Bjjfc given from Eqs. ( p9\) and 
( [zol ) but with V{r) replaced by T2b(eo^ + e^^) in all matrix elements. AE^{p) corresponds to the perturbative, 
two-body limit of Ai?3 (p) and the expression for AE!^ {p) is finite even if the contact potential approximation is used 
for T2b. 

The physical interpretation of the renormalized quantities and AE^{p) is that they describe the difference 
between scattering in a medium and scattering in a vacuum. They therefore upgrade the two-body T-matrix which 
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now appears in leading order matrix elements and replace it with the many-body T-matrix. The use of ordinary 
perturbation theory on these terms corresponds to a second order Born approximation of Tmb in terms of [c.f. 
Eq. (pO|)]. In contrast to a Born expansion in ^(r), this is expected to be valid for a low temperature, dilute gas away 
from the critical point. As before, self-consistent perturbation theory introduces Tmb to all orders in condensate- 
condensate collisions but still only to second order in condensate-excited state collisions. Since the expansion is now 
in terms of T2b, however, the difference between these is of higher order than the calculation of this paper. Thus, now 
that matrix elements arc written in terms of the two-body T-matrix, a perturbative treatment of the non-quadratic 
terms in the Hamiltonian should be justified, at least away from the critical point. The validity of the theory in this 



form is discussed further in Sec. VIl 



Although we have given explicit formulae for the renormalization of and Ai^a, it is often sufficient simply to 
neglect the zero-temperature contribution to these quantities. This part of the expressions contains all the ultra- 
violet divergence because the quasiparticle populations which appear at finite temperature vanish exponentially at 
high energy. Of course this procedure is not exact as it neglects the difference between the many-body and two-body 
T-matrices at zero temperature. For dilute gases these corrections are of order (na"^)^/^ and so can justifiably be 
ignored for the trapped gases of current experiments. 



2. The ground state energy 



In the previous section we have shown that the use of the two-body T-matrix to describe particle interactions 
leads to finite expressions in the GPE and BdG equations. It remains to show that the total energy of the system 
(and in particular the ground state energy) is also finite as this will establish that the theory is completely free of 
ultra-violet divergences. The renormalization is rather more subtle in this case, however, and the contribution from 
the non-quadratic Hamiltonian is logarithmically divergent if the contact potential is used. 

The total energy of the system is given to the order of this calculation by 



E — Hq + E2 + E^ 



E4 



(109) 



Writing this in terms of the two-body T-matrix and introducing the appropriate ultra-violet renormalization gives 
[c.f. Eqs. d), (H, (H) and (||)] 



E = No[h, 



(00|T2b|00) 



i^.T + 27Vo(0^|T2b|J0)|p,, 



^ i^(oo|T2bKj)4 + ^(.j|r2b|oo)4 



ij^O kin 



(0»|T2b|fcm)^P^P(fcm|r2bb-0) 



(110) 



4 51 («il^2b|fcm) [ 



PkiPmj + PmiPkj + '^ij* l-^mk] 



where for simplicity we have neglected the z dependence of the T-matrix. We note that only one of the linear 
contributions from in this equation is renormalized and the other still contains an ultra-violet divergence if we use 
a contact potential. The single-particle energy X^y/o ^^^^ divergent, however, and it turns out that the two 

divergences cancel exactly. 

The ground state energy is obtained by using the zero-temperature expressions for ptj, Kij and E3. The renormal- 
ization leads to a finite contribution at quadratic order and in the homogeneous limit we obtain the well-known result 



E, 



a 
N 



2 



128 
150F 



(111) 



where n is the total number density. At next order we need to consider the effect of E^ and -E4 . The renormalization 
of Ei leads to a finite contribution of order na^ , while the renormalization of E^ removes its linear divergence, but 
leaves a logarithmic term. This can be dealt with using the fact that the full two-body T-matrix has a high-energy 
cut-off around fc ~ 1/a. Truncating the summations for E^ in this region leads to an expression for the ground state 
energy which has been calculated by Hugenholtz and Pines |l| and Wu PI 
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a 
N 



2 



1 + J^(^a3)i/2 + 8 (l^i „ y3)(na3) in(na3) + 0{na^) 

ISVTT 3 



(112) 



The term of order no?' in this expression depends on the exact position of the cut-off and hence on the particular form 
of the interatomic potential. Thus to the order of the calculation we have presented, the contribution to the ground 
state energy at order na^ is the only quantity which can not be calculated solely in terms of the s-wave scattering 
length. 



VI. INFRA-RED DIVERGENCES AND A GAPLESS SPECTRUM 



As mentioned in the introduction, two difficulties commonly associated with attempts to extend the Bogoliubov 
method beyond the quadratic Hamiltonian are the presence of a gap in the excitation spectrum at low energy and the 
appearance of infra-red divergences in the higher order energy shifts. Genuine infinities do not appear in a trapped 
gas because this has a natural low energy cut-off set by the trap frequency. Divergences may re-appear, however, as 
the trap is opened and the system approaches the homogeneous limit. Thus for a trapped gas the problem of infra-red 
divergences corresponds to an anomalous dependence of physical quantities on the trap frequency. A discussion of 
infra-red divergences and a gapless spectrum is therefore a discussion of the homogeneous limit, and if the theory is 
well-defined in this case we can assume that this is also true in a trap. 

In this section we will show that the theory developed in this paper is finite and gapless in the homogeneous limit. 
We will restrict our attention to neutral atoms in three dimensions for which the contact potential can be used as an 
approximation to the two-body T-matrix. Details of the relevant calculations can be found in Ref . . Our notation 
follows that of Refs. |l|,|^. 



A. Quadratic theory 



In the homogeneous limit, the GPE of Eq. Q [with V{v) T2h{z) Uq S{r)] leads to a condensate wave function 
and eigenvalue given by Co = 1/v^ and A = uqUo where no = Nq/H and Q is the volume of the system. The wave 
functions Ci(r) which describe the noncondensate can be chosen to be plane waves Ci(r) = since these 

form a complete set orthogonal to the condensate. The coefficients Cij and Mij of Eq. ( p^ ) which define the BdG 
equations are therefore 



Cij — 



2 1,2 



2m 



uoUq 



where ki = |ki| and the label —j refers to the state with wave vector k_j 
The quasiparticle transformation of Eq. (|l^) reduces to 



Pk = Ukcik - Vka^_k, 



and the solution of the BdG equations gives 



1 



Uk 



where 

ak = l + 

and the dimensionless wave vector yk is defined by 



Vk 



-Oik 



Vk 



Vk = ^, 
fco 



2m 



= noUo ^ fcg = SnauQ. 



The quasiparticle energies are given by the usual Bogoliubov expression 



2m 



1/2 



noUo yk{2 + ylY^'^ = uqUq Zk, 



(113) 



(114) 



(115) 



(116) 



(117) 



(118) 
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where c = {noUo/mY^'^ is the speed of sound, and the final term on the right hand side simply serves to define the 
dimensionless energy Zk- 



The expression for the quasiparticle energies of Eq. (118) must be supplemented by an expression for the condensate 
density in terms of the total number of particles. This can be obtained by solving the nonlinear equation n = uq — p^x, 
where the noncondensate density pex is obtained from the expression [c.f. Eq. (44)] 



1.3 

(27r)= 



d^Vk 



nk + {n^k + 1) al 



I -at 



(119) 



For a dilute gas at zero temperature, the condensate population is therefore given to order (na^Y^^ by the well-known 
result in 



no(T = 0) ~ n 



3\ 1/2- 



The renormalized anomalous average of Eq. (103) is given by the expression [c.f. Eq. (45)] 



1,3 

" (27r)3 



{uk + + 1) Ok 1_ 

2yl 



(120) 



(121) 



where the last term corresponds precisely to the ultra-violet renormalization of Eq. (|10J). n"' is finite at any temper- 
ature and at T = we have the result k" — 3pex — fcg/(2V27r^). 



B. Higher Order Terms 



For calculations beyond quadratic order, the properties of the condensate are calculated from the generalized GPE 
of Eq. (IOC) rather than the ordinary GPE. The condensate wave function is unchanged, but the eigenvalue becomes 



Ag = n-oUo + 2UoPex + Uqk'^ 
Evaluating this at T = 0, we obtain the well-known result jsj 
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AG(r = 0) = n;7o 



1 



(122) 



(123) 



The discussion of Sec. IV shows that the total contribution to the quasiparticle energy from non-quadratic terms 
can be written as 



A£:(fc) = Ai;3(fc) + AEiik) + AEx{k). 



(124) 



The energy shifts Ai?4(fc) and AEx{k) are defined in Eqs. ( |6l| ) and (|64|). Using these, together with Eq. ( |122| ), we 
obtain 



(125) 



This result demonstrates that in the homogeneous limit there is no shift in the quasiparticle excitation spectrum in 
the HFB-Popov approximation, where both Ai?3(fc) and are neglected. This result is expected from our discussion 
of Sec. p^VE where we showed that the energy shifts in the Popov approximation depend on the variation of the 
noncondensate density in the region of the condensate. 

The energy shift AE^^k) is defined in Eq. (74) and can be written in the homogeneous limit as 



AE3{k) = - 



Uokl 



[AE^k] + AE',{k) + AElik)] , 



where 



AE-{k) = / d^y, 



{1 + rii + n'j).[l — Qfj — aj + a^a^ + oijOik ~ aiOijOn^ 



(126) 



(127) 
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(1 



(z, + Zj- + Zfc).(l-a|).(l-a2) 



ivl - yf - yj) 



(128) 



AE-^ik) = -2 / d^y, 



(129) 



(i = k - j) 



The final terms of Eqs. (127) and ( |128| ) correspond to the ultra-violet renormalization of Ai?3(fc) given in Eq. (lOS). In 
the above equations a mixed notation of a's, z's and j/'s has been used in order to keep the expressions as compact as 
possible. The evaluation of the integrals is greatly simplified using the following relationships between these quantities 



211/2 



ttfc = [1 + 

1 - "fe = 2zfeafc, 



Zk, 



211/2 



k' = [^ + 4] 

(1 - akf = "^yiak- 



Zk, 



(130) 



The integrals of Eqs. (127)-(129) contain energy denominators which vanish when a decay process is energetically 
allowed. The integrals must therefore be calculated by inserting a small imaginary part le and using the result 



lim — - — = V {- 

e— »o X + le \x 



nrS{x), 



(131) 



where V means 'the principal part'. The first contribution describes the energy shifts of the excitations, while the 
imaginary part corresponds to their decay and describes the quasiparticle lifetimes. These two contributions are 
discussed in the following subsections. 



1. Energy shifts 



The energy shifts at low momentum can be calculated by introducing an expansion in powers of zj, [we could equally 
well use an expansion in powers of since for <C 1 we have Zk — \/2yk + 0{y^)]. This expansion has the form 



AE: 



3, 4, A 



(fc) 



A 

Zk 



B + Czk + 0(zl 



(132) 



The origi n of the te rm in 1/zfe is the prefactor 1/(1 — a|), which appears in each of AE-i{k), AE4{k) and AE\{k) 
[c.f. Eqs. (125) and (126)] and which has the expansion 



1 

2zk 



Zk 



+ 0{zl). 



(133) 



The presence of terms proportional to l/zk means that there is an infra-red divergence in both A£^3(fc) and AE4{k) + 
AE\{k). To prove that there are no infra-red divergences in the overall energy shift, we therefore need to show that 
when we add these quantities the resulting coefhcient A is zero. To show in addition that the theory is gapless we 
need to show that the sum of the B coefficients is also zero. The total C coefficient then describes the energy shift of 
the low-momentum states and corresponds to a modification of the speed of sound. This should be small if the use of 

perturbation theory on the non-quadratic Hamiltonian is justified. 

The expansion of AEi{k) + AE\{k) can be obtained directly from Eq. (125) and is 



AEi{k) ^ AEx{,k) ^ ~ 



2C/ok" Uon'^Zk 
Zk 2 



0(zD. 



(134) 



This result demonstrates the presence of infra-red divergences in the full HFB theory (where AE^{k) is neglected) 
which we mentioned earlier. The calculation of AE^{k) is more difficult and the details are given in Ref. ||35| . The 
leading order contribution is simply obtained by setting i = j in Eqs. (127)-(|l2^). This gives 



AE^,{k) = AElik) 



hm -TT-^r^ — 



d^y, 



1 



2zi 



2tiY 



(135) 
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while AE^{k) = 0. Comparison with Eq. (121) [using the fact that a/c/(l — a1) = l/{2zk) from Eq. ( |l30| )] gives 



AEsik) = AElik) + AE^ik) = --^ + . . . (136) 

Zk 



Equations ( 134 ) and (|136| ) show that the infra-red divergence in AE3{k) cancels exactly with that in Ai?4(fc) + 
AEx{k). We note that this occurs without the need to assume any particular functional form for the quasiparticle 
populations. Thus the infra-red divergences cancel for non-equilibrium quasiparticle distributions as well as for thermal 
distributions of any temperature. This general result is what is needed to show that the theory is well-defined, since 
the cancellation of divergences should be an intrinsic property of the theory and not of a particular distribution. 

To show in addition that the theory is gapless requires a calculation of A£'3(fc) to order (zfc)". T he re sult is that 



this term vanishes, and since there is also no contribution at this order from AE4^{k) + AE\{k) of Eq. (134) this means 
that the theory is gapless. We have proved this result for any quasiparticle distribution which depends only on the 
modulus of the wave vector and which varies at worst as Uk ^ ^/Vk as yk — > 0. Since the low-momentum excitations 
are phonons with a linear dispersion relation, this includes the Planck distribution at finite temperature. The overall 
shift to the quasiparticle energies from the non-quadratic Hamiltonian is of order Zk, and the theory is therefore finite 
and gapless. 

Thus the effect of the non-quadratic terms on the phonon spectrum is simply to modify the speed of sound. At 
zero temperature we have calculated the energy shift to be 

AE{k, T = 0) = Zk + 0{zl). (137) 

This shift is small compared to the leading order expression = noUoZk (justifying the use of perturbation theory) 
if the dilute gas criterion na^ ^ 1 is satisfied. At high temperature {kBT/noUo ^ 1), the energy shift has been 
calculated by Fedichev and Shlyapnikov |jri| and is 

AE{k, T) « -In^UoZk {n^a^f'\ (138) 



which is small if [ksT /ntiUQ).{nQaP'Y^'^ 1. This i s the refore the condition for the validity of the theory in the high 



temperature regime and is discussed further in Sec. VII 



2. Decay rates 

AE3{k) develops an imaginary part AEl{k) whenever a real decay process is energetically allowed. T he c orre- 
s pon ding decay rate 7(fc) = —2AE\{k) /Ti can be calculated by replacing the en ergy denominators in Eqs. (127) and 
( |129| ) with ~tn6{E), where S{E) is an energy conserving delta function [c.f. Eq. (131)]. This leads to the same results 
as the application of Fermi's Golden Rule to the cubic Hamiltonian for the case that the decaying mode has a large 
population. 

Integrating over the energy conserving delta function, we obtain an expression for the decay rate in the form 
j{k) = 7s (fc) -t- jL{k), where the Beliaev decay rate jsik) is given by 

hknU f^'' (1 + Hi + n^).\l — Q!j — tto -|- Ckiaj. + a^ak — aia^akV 
^B[k) — / dzj 



SmykZkakJo ' a,aj{l + zfy/^{l + z^y/^ 

(zk ^ Zi + Zj), (139) 



and the Landau decay rate 7L(fc) is given by 



hk^a f°° [n-j - ni).[l - Uj - ak + aiUk + aiUj - aiUjakY 



^^^''^ SmykZkakJo '^'^ a,a,il + zfy/^1 + zjy/^ 

{zk = z,-Zj). (140) 

Beliaev damping corresponds to processes in which the decaying quasiparticle decomposes into two of lower energy, 
while Landau damping is essentially the reverse process in which two quasiparticles collide and annihilate to form a 
third of higher energy. 
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Using Eqs. ( |139| ) and ( |140| ) all the standard results for the decay rates of a homogeneous Bose gas can be obtained. 
For example, at T = the decay rate is entirely due to Beliaev processes and for the phonon part of the spectrum it 
is given by [p 17|li 



At high temperatures Landau processes dominate and the phonon decay rate is given by ||ll|.y_2|,^ljj42 



7 ... 

2 h \noUo 



^ {noa^)^/\ {kBT/noUo » 1, fc/fco « !)■ (142) 



C. Physical Discussion 

We have shown that the inclusion of A£'3(fc) leads to the cancellation of infra-red divergences and to a gapless 
spectrum. We will now discuss the physical reasons why these difhculties can arise in the theory of EEC and why 
individual terms in the perturbation theory contain divergences. 

The excitation spectrum is calculated as the energy required to create a quasiparticle. This involves a change in the 
number of excited atoms which is A{Nex) = (1 + «fc)/(l — c^l) ^ hmfe^o l/(V22/fc)- ^ single low-energy quasiparticle 
therefore contains a very large number of particles and this is the physical origin of infra-red divergences in individual 
terms of the perturbation theory. A consequence of this is that care must be taken to ensure that the total number 
of particles is fixed so that the additional excited particles are taken from the condensate. This is the origin of the 
contribution from AE\{k) which is crucial in producing a well-defined theory. 

The fact that low-energy quasiparticles consist of a large number of particles can also be seen as the reason why 
the low- momentum excitations of the system have the phonon dispersion relation = chk. This involves a non- 
perturbative change in the excitation spectrum from the quadratic form appropriate to particles, and occurs because 
the energy of a quasiparticle corresponds roughly to the energy of a particle multiplied by the number of particles it 
contains. Since the latter quantity scales as 1/fc as fc this gives Cfc ~ (?i^fc^/2m).(fco/\/2fc) = chk/2, which only 
differs from the actual value by a factor of two. 

The mathematical origin of infra-red divergences is the fact that the quasiparticle functions Uk and Vk scale as 1/fc 
for fc — > 0. The energy of a quasiparticle is a quadratic form in the u's and ?;'s and the only quadratic combination 
which is finite as fc ^ is — = 1. Solving the BdG equations ensures that the quasiparticle energies depend 
only on this particular combination and hence are well-defined. If we modify these equations, however, and treat the 
changes via ordinary perturbation theory, then we will get energy shifts which are still quadratic forms in the it's 



and v's but which are not necessarily proportional to — [c.f. Eq. (125)]. In this case there may be an infra-red 
divergence simply because of the low-momentum behaviour of the u's and w's and not because of the particular change 
to the BdG equations. Explicit calculations as in the previous section are required to show whether or not infra-red 
divergences cancel in physical quantities. 

If we use self-consistent perturbation theory, however, and diagonalize the BdG equations exactly then these infra- 
red divergences will disappear, although the remnant of the problem will be seen in the energy spectrum having a gap 
at low energy (as in the HFB theory). Self-consistency therefore leads to finite energy shifts, but these will certainly 
be large compared with the quadratic theory results if the approach is not gapless. Any such theory is therefore 
inconsistent because a perturbative treatment of the non-quadratic Hamiltonian requires that the energy shifts are 
small compared to the leading order terms. 

The appearance of infra-red divergences (and hence a gap in the excitation spectrum) can also be seen as a 
consequence of an inconsistent treatment of interactions involving the condensate. If Ai?3(fc) is neglected (as it is in 
the HFB theory), then the many-body T-matrix is introduced into pair excitation processes of the form |00) 
but not into the Hartree-Fock interactions (Oz|y''|jO). These terms are of the same order, however, as each involves 
two condensate labels and two excited state labels. The different treatment of these interactions in the HFB theory is 
therefore inconsistent and is responsible for the gap in the excitation spectrum. In contrast, the inclusion of AEt, leads 
to a consistent theory and a gapless spectrum. We showed in Sec. 0, however, that AE^ only introduces the T-matrix 
to second order in the Hartree-Fock terms, whereas self-consistent perturbation theory with the anomalous average 
gives the T-matrix to all orders in pair excitation processes. This means that the theory developed in this paper may 
not be gapless if self-consistent perturbation theory is used. Since we have proved that ordinary perturbation theory 
gives a gapless spectrum, however, the size of any gap will be beyond the order of this calculation. In contrast, the 
size of the gap in the HFB theory is proportional to and hence is of the same order of magnitude as the calculation 
which is being performed. 
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The above argument also shows why the Popov approximation produces a gapless theory. There are essentiahy 
two consistent ways to treat condensate interactions; either we include the effects of the medium in all collisions or 
we neglect them completely. The first approach corresponds to the theory discussed in this paper and requires the 
inclusion of lS.E^[k). The second approach corresponds to the Popov approximation in which both Ai?3(fc) and n'^ are 
neglected. In this case all interactions are desc ribe d by the two-body T-matrix (contact potential). The gapless nature 
of the Popov approximation is shown by Eq. ( [l34|) which also demonstrates that there are no infra-red divergences in 
this approach. Thus the HFB-Popov the ory is p referable to the full HFB approach and for this reason it has been used 
in recent calculations for trapped gases ||3y,^. Nonetheless, the neglect of many-body effects can only be justified 
at near zero temperature and the full second order theory of this paper will be needed whenever there is significant 
depletion of the condensate. 



VII. VALIDITY 



The principle assumption we have made in this paper is that the existence of a condensate means that the non- 
quadratic terms in the Hamiltonian can be treated using perturbation theory. We argued in Sec. 0, however, that 
this can only be the case if interactions are described by the two-body T-matrix rather than the interaction potential. 
The discussion in this section will assume that this has already been done. 

The use of perturbation theory requires that the effect of non-quadratic terms is small, and this can be verified after 
any calculation. A necessary condition is that the total energy shifts do not contain infra-red divergences and vanish 



at low momentum. We have shown in Sec. VI that this is indeed the case and the validity of the theory therefore 



depe nds on the size of the resulting shifts relative to the leading order terms. A comparison of Eqs. (US), (137) and 



(138) shows that the validity of the theory in the homogeneous limit can be summarized by 

(na3)i/2<l, (T = 0), (143) 



(noay/^ « 1, {ksT/noUQ » 1). (144) 



These same conditions are obtained by requiring that the quasipar ticle s are well- defined (i.e. that their decay rates 



are small compared to th eir f requencies) as can be seen from Eqs. ( |141| ) and (142). 

The condition of Eq. (|l43| ) for the validity of the theory at zero temperature is simply the criterion for a dilute 
gas. The finite temperature condition of Eq. ( |144| ) is closely related to the Ginzburg criterion for the validity of 
mean-field theory [Q. The fact that this expression scales with the condensate density as (rio)~^^^ means that the 
theory must fail in the region of the critical point where — > 0. The boundary of the critical region can be estimated 
by setting {kBT/noUa){noa'^y^^ ^ 1. If we write this in terms of the dc Broglie wave length of the atoms and use the 
noninteracting gas result for the condensate density near the critical point, then it becomes t^/^ ~ a"'^/^/(n^/^A^^), 
where t is the reduced temperature defined by t = (T^ — T)/Tc for T < Tf.. Since nX^g ~ 1 near the critical point, this 
can be simplified to t ^ a/XdB- The critical region is therefore very small for a dilute gas for which a/Xds ^ 1- We 
note that the failure of the theory at the critical point is a result of the growth of long wavelength fluctuations which 
destroy the mean field of the condensate. This 'infra-red' behaviour produces a genuine breakdown in perturbation 
theory and leads to the occurrence of critical phenomena. It is completely distinct, however, from the problem of 
infra-red divergences in individual terms of the perturbative expressions with which we have been concerned in this 
paper (and which we ha ve sh own cancel in physical quantities) . 



Equations (143) and (144) show that the expansion parameter in the homogeneous, thermodynamic limit is pro- 



portional to a'^/'^ at T — and a}/'^ at finite temperature. These are non-analytic functions of a which means that 
the perturbation expansion is at best only asymptotically convergent. This is demonstrated both by the presence of a 



logarithmic term in the expression for the ground state energy of Eq. (112), and by the fact that a condensate can not 
exist with a negative scattering length |44| . In a trap, however, the finite size of the system means that a condensate 
can exist in the negative scattering length case, provided that its population is not too large [[l5|j46| ] . This leads to 
the possibility that the expansion is convergent in a trap, at least for small condensate populations. 

This result indicates that the trap potential can have a profound effect on the nature of BEG and the validity of 
mean- field theory. In particular, it is possible that the theory developed in this paper can be meaningfully applied 
to a trapped gas even in the region of the phase transition, at least for some range of values of the s-wave scattering 
length. The reason is that the finite length scale in a trap suppresses the long wavelength fluctuations which are 
responsible for destroying the condensate mean held. One obvious consequence of this is that BEG is possible in one 
and two dimensions in a trap, whereas in the homogeneous, thermodynamic limit it can not occur in fewer than three 
dimensions. 
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A criterion for the validity of the theory which can be apphed to a trapped gas is provided by the anomalous average 
Kij . This is an important quantity because it vanishes in the absence of interactions and is only large when there is a 
significant noncondensate fraction. The anomalous average introduces many-body effects into particle collisions and 
the requirement that these are small compared to two-body effects can be written as 

E.,,o(00|^.N.)^j^ « 1. (145) 



In the homogeneous limit this becomes K^/no <C 1 and leads directly to the conditions of Eqs. (143) and (144). The 
fact that the finite temperature criterion can not be satisfied close to the phase transition can be seen from the fact 
that the many-body T-matrix vanishes at the critical point for a homogeneous gas ||l^,^. In this case the effects 
of the medium on the collision of two particles are as large as the two-body effects and the pcrturbative expansion 
breaks down. In a trap, howeve r, the many-body T-matrix remains finite in the region of the phase transition and 



only decreases by of order 10% |33 4^. This is a further indication that the theory may remain valid through the 
phase transition in a trap. 

Another simple test of the theory in a harmonic trap is provided by the calculation of the Kohn modes [Q . These 
are excitations which consist of a centre of mass oscillation of all the atoms about the trap minimum, and they are 
exact solutions of the equations of motion for an interacting system. We therefore expect there to be solutions of the 
BdG equations with frequencies of exactly one in the trap units appropriate to the axis of oscillation. 

The Kohn modes are recovered exactly in the quadratic theory because this includes the full dynamics of the 
condensate. They are not obtained exactly in higher order theories, however, because the motion of the thermal cloud 
is treated approximately. This can be seen from the fact that first order perturbation theory is used on the quartic 
Hamiltonian _ff4, corresponding to the assumption that the noncondensate is static. Some of the dynamics of the 
thermal cloud are included in this paper, however, because second order perturbation theory is used on the cubic 
Hamiltonian H^. Nonetheless we do not expect to recover the Kohn modes exactly. This fact can be turned into a 
useful estimate of the error in the theory, however, since our method is based on a systematic treatment of the original 
Hamiltonian, so the degree to which an exact result is violated provides a measure of the importance of higher order 
terms. 

Finally, we mention that although our discussion of validity has focused on the regime below the critical temperature, 
the leading order predictions of the theory above Tc should also be correct. The reason is that in the zero-condensate 
limit the theory reduces to the Hartree-Fock approximation for the noncondensate, which should be reasonably 
accurate for a dilute gas above the critical point. In addition, the leading order properties of the system at high 
temperature are simply those of a noninteracting gas, and this contribution is treated exactly in the quadratic 
Hamiltonian. For this reason numerical simulations for trapped gases reproduce the noninteracting gas results for 
r>Tc [|0|j3l]j3|]. 



VIII. GAPLESS HFB 

In this section we describe and motivate a gapless extension of the conventional HFB theory which has recently 
been developed [^,0. This theory has the advantage that it is significantly easier to compute with than the full 
second order theory of this paper, and it can therefore be used as a first approximation when many-body effects are 
not negligible (so that the HFB-Popov approximation can not be applied). 



In Sec. VI we argued that the physical reason for the appearance of a gap in the HFB theory is the inconsistent 
treatment of interactions involving the condensate. In particular, the anomalous average introduces the many-body 
T-matrix into the off-diagonal terms of the BdG equations, but if Ai?3 is neglected then the diagonal terms only 
contain the two-body T-matrix. A simple way to solve this problem is therefore to insert the many-body T-matrix 
'by hand' into the diagonal terms by replacing the matrix elements (Oi|T2b|iO) with (Oi|Tmb|iO). This is the basis of 
the gapless HFB theory and it is also equivalent to the effective Hamiltonian approach of Bijlsma and Stoof p^ . 

If the two-body T-matrix is approximated by a contact potential then the low-energy limit of the many-body 
T-matrix will also have this form. The diagonal matrix elements can therefore be written as 



(0*|r„,bb-0) = jd\C*{v)Q{v)T^.y.{v)\Co{r)\\ (146) 
where we have denoted the T-matrix by Tmblr), since many-body effects are spatially dependent for a trapped gas. 
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The interaction strength Tinb(r) can be found from consideration of the off-dia gon al matrix elements^ («i|Ti„b|00). 
Writing Eq. ( ^o| ) in the position representation and using Eqs. (p3|), (101) and (|l02| ) we have0 



No /d3rC(r)C;(r)T„,b(r)Co'(r)= /d^r C (r)C; (r) [iVoC/oCo W + C/o^^'W] 



(147) 



No /d^rC*(r)C;(r)[/o 



K«(r) 



A^oCo'(r) 



We can therefore obtain a gapless HFB theory by starting with the generalized HFB approach (i.e. K"(r) ^ 0, 
AE3 = 0) and introducing the many-body T-matrix to all orders in the diagonal elements of the BdG equations by 
replacing Uo with Tnib(r), defined by 



r,nb(r) = Uo 



K«(r) 



(148) 



The result of Eq. (^9|) shows that Tmh (r) is real and therefore the quasiparticles do not decay in this approach. 

The discussion to this point has focused on the condensate-excited state collisions (OilTmbljO) which are the leading 
order terms in the diagonal elements of the BdG equations. The introduction of the many-body T-matrix in these 
interactions is mandatory if we wish to have a gapless theory. However, there are further interaction terms in the 
BdG and GPE equations which only contain the two-body T-matrix at the order of this calculation. In particular, 
collisions of noncondensate atoms are described in the BdG equations by the term X^/cm^^o ^(^*l'^2b|j'7i)pm/c which 
appears in Eq. ( |l05| ), while the effect of condensate-excited state collisions o n th e condensate is described by the 
expression X^ij^^o '^{ki\T2h\jO)pji which appears in the generalized GPE of Eq. ( |100| ).p| The theory remains gapless if 
these interactions are also upgraded to the many-body T-matrix and we therefore obtain two versions of the theory 
depending on the treatment of these terms. In the first approach (GHFBl) these collisions remain described by the 
two-body T-matrix [/q, while in the second approach (GIIFB2) many-body effects are introduced by replacing Uo 
with rinb(r). 

The two gapless HFB theories can therefore be summarized by the following equationsFj 



'Trap(r) 



^ Afor„,b(r)|Co(r)PCo(r) +2C/e,(r)pe.(r)Co(r) = XcCoir), 



(149) 



£{v)u,{v) + NoTn,^{r)Cf>{r)v,{v) = e,u,(r), 

(150) 

£(r)i;,(r)+iVorn,b(r)Co*'(r)w,(r) = ~e,v,{r), 



where £(r) is defined by 



h 

C{r) = -— + VTrap(r) -Xg + 2A^oT„,b(r)|Co(r)P + 2U,^{r)p,,{r), (151) 
zm 

and we have used the fact that rmb(r) is real. The GHFBl theory corresponds to setting C/cx(r) = '^o while the 
GHFB2 theory is defined by J7cx(r) — 7inb(r)- We note that there may be numerical difhculties in the implementation 
of the GHFB2 theory because the BdG equations contain the term 2Uo [l+K"-{r)/No (oi'"^)] Pea;(r), which may not 
be well-defined at the edge of the condensate. There is no such problem in the GHFBl theory because in that case 
Tnib(r) always appears multiplied by Co(i")- 



'^^We are assuming that the dependence of the T-matrix on the energy of the collision (the parameter z) can be neglected 
at low-energy. In reality the many-body T-matrix does have a momentum dependence and at high-energy it reduces to the 
two-body T-matrix. This ultimately leads to two versions of the gap less HFB theory as discussed later in this section. 

^''We note that is defined here using Tj^b rather than T2I) in Eq. (102). This definition is more appropriate to self-consistent 
perturbation theory, but the two are equivalent to the order of this paper. 

^""The interaction strength which appears in terms involving the anomalous average must remain written in te rms of r2b, 
since this is what is required to introduce the many-body T-matrix into condensate-condensate collisions [c.f. Eq. ( |l47|)] . 

^®To obtain excitations orthogonal to the condensate we should really write these equations in the form of Eq. (|38[) with 
Uo rmb(i') in Et- (0)- Since the theory is gapless, however, we can simply apply the method of Eq. (^^. 
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The additional many-body effects introduced in the GHFB2 theory can be justified on the physical grounds that 
all low-energy collisions should be described by the same effective interaction. Nonetheless, the many-body T-matrix 
depends on the relative momentum of the colliding particles and Ti„b(r) is only a good approximation at low energy. 
The GHFB2 theory is therefore only appropriate if most of the collisions involving the noncondensate occur at low 
energy. If the dominant collisions occur at an energy which is large compared to the mean-field interaction, then 
the GHFBl theory should be used because in this limit the many-body T-matrix reduces to the two-body T-matrix. 
The two theories therefore correspond to the two extreme assumptions about the energy regime of noncondensate 
collisions. 

We will now discuss the relationship between the gapless HFB approach and the full second order theory of this 
paper. The essential feature of the GHFB theories is the introduction of the many-body T-matrix into condensate- 
excited state collisions in the diagonal terms of the BdG equations. We showed in Sec. ^that such T-matrix corrections 
are introduced by albeit only to second order. AEt, also introduces polarization effects, however, and these are 

not contained in the many-body T-matrix. The GHFBl theory can therefore be seen as an approximation to the full 
theory developed in this paper which neglects these polarization terms, but includes higher order contributions to the 
T-matrix. The additional many-body effects in the GHFB2 theory are also of higher order and will be introduced in 
a full theory by higher order perturbation theory as discussed in Ref. ]35| . The higher order terms in both the GHFB 
theories should not have a significant effect, however, and if they are important then a systematic calculation beyond 
the order of this paper will be required. The gapless HFB approach is only useful therefore if the dominant contribution 
from AE^ corresponds to many-body T-matrix corrections, and if terms beyond the order of this calculation can be 
neglected. 

The quasiparticle energies predicted by the gapless HFB theories can easily be calculated in the homogeneous limit. 
In this case, the total energy shift from non-quadratic terms can be written as 



AEik) = A£;ghfb(A:) + AE^ik) + AExik), 



(152) 



where Ai?4(/c) -f- AEx{k) is given in Eq. (125) and the contribution from the terms introduced in the GHFB theories 
is 



AE, 



GHFB 



(fc) = 2J7ok" 



(153) 



with ak and k'^ defined in Eqs. (116) and (121) respectively. We note that the two gapless HFB theories give exactly 



the same energy shifts in the homogeneous limit because the additional many-body effects introduced in GHFB2 
cancel with each other. 



Substituting Eqs. (125) and (053) into Eq. (152) gives 



AE{k) = Uqk 



In the low-momentum limit fc ^ this becomes 



AE{k) 



0{zl) 



(154) 



(155) 



Evaluating using Eq. (121), we find that the zero and finite-temperature energy shifts for the phonon spectrum 
(fc -C fco) are 



AE{k,T = 0) = 



Zk, 



(156) 



AE{k,T) = -2tt 



'^^noUoZk (^) (noa^)l/^ (fcsT/not/o » 1). 

\noUoJ 



(157) 



Comparison with Eqs. (137) and ( 138 ) shows that these results agree with those from the full theory to within factors 
of order unity. 

For trapped gases, the gapless HFB theories have been implemented numerically in both one and three dimensions 
[^,|^. In one dimension, the predictions of GHFBl for the quasiparticle energies did not differ significantly from 
those of the HFB-Popov theory, although there was a noticeable difference in GHFB2 . This is somewhat surprising 
in the light of the homogeneous limit results (where all the effects occur in the GHFBl theory), and the fact that the 
additional many-body corrections introduced in GHFB2 are expected to be of higher order. In three dimensions, the 
GHFB2 theory gives good agreement with the experimental results for one of the low- lying excitations at temperatures 
approaching Tc, although the agreement for another is worse than the HFB-Popov predictions |p4]. 
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IX. CONCLUSIONS 



In this paper we have developed a gapless theory of BEC which goes beyond the quadratic theory of Bogohubov 
in a consistent manner, and which can be apphed to trapped gases as well as to the homogeneous limit. The theory 
is based on rewriting the many-body Hamiltonian in a form which is approximately quadratic, and using first and 
second order perturbation theory to treat the non-quadratic terms. Infra-divergences can appear in individual terms 
of the perturbation expansion, but the total contribution is finite and the change in the excitation spectrum is small 
for a dilute gas away from the critical point, justifying the use of perturbation theory. The problem of ultra-violet 
divergences is dealt with by using the contact potential as an approximation to the two-body T-matrix rather than 
the bare interaction potential. This leads naturally to a renormalization of the theory at high energy. 

The numerical implementation of the theory for trapped gases is currently in progress |39[| , and it is hoped that this 
will lead to good agreement with experimental results for the energies and lifetimes of the excitations at temperatures 
approaching the phase transition. The theory can also be applied to the study of charged Bose gases and neutral 



atoms in two dimensions, and work along these lines is also in progress 49 1. 

A further are a of interest for future research is the study of the properties of BEC near the phase transition. We 
showed in Sec. VII that the theory developed in this paper fails in the homogeneous limit near the critical point, 
but may remain valid in this region for a trapped gas. Whether or not this is the case can be established using the 
requirement that the perturbative shifts predicted by the theory are small. The signature for the failure of the theory 
is the appearance of large shifts for low-energy states, and this should occur as the trap is opened and the system 
behaves in a quasi-homogeneous manner. 

The failure of perturbation theory is associated with the emergence of critical phenomena, and is a consequence of 
the growth of long wavelength fluctuations on the condensate. This can only be studied using an approach which is 
free of the spurious infra-red divergences which arise from an inconsistent treatment of the non-quadratic Hamiltonian. 
It is therefore important to apply the theory developed in this paper to the region of the phase transition and to 
study the cross-over between mean-field and critical behaviour for a Bose gas. 
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APPENDIX A: EFFECTIVE LOW-ENERGY HAMILTONIAN 



In this appendix we will show for a Bose gas at low temperature that the high-energy states of the system can be 
adiabatically eliminated to give an effective Hamiltonian for the low-energy states. The interactions between particles 
in the low-energy subspace are described by an approximation to the two-body T-matrix. 

We start from the many-body Hamiltonian for a system of structureless bosons with pairwise interactions [c.f. 
Eq. &)] 



H 



E 

ijkm 



{ij I V I km) dj OjCikO,. 



(Al) 



huJiSij and the index i refers to 



where for convenience we have used the basis in which Hgp is diagonal,]^ so H, 

the wave function (^^^(r). The operators aj and Ui are respectively the creation and annihilation operators for the 
eigenstates of H'^^ and have the usual Bose commutation relations 



0. 



(A2) 



We will subsequently make a unitary transformation to the basis orthogonal to the condensate which is used in the main 
text. 
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The equation of motion for any operator O is given by the Heisenberg formula 



in—r- 
dt 



0,H 



Using this we obtain the equation of motion for the annihilation operator 

. da 
~dt 



jkm 



(A3) 



(A4) 



where uj. 



ijkm 



uji + (jjj — ujk — ujm and we have eliminated the free evolution by defining at 



on products of three operators so we will also need the equation of motion for these, which is 



Eq. (A4) depends 



d{a\akam) v-^ -t- 
ih — = 2, (fcm|y |pq)ajapaqe" 



dt 



(A5) 



pq 



pqr 



{pm\V'\qr)a]alakaqare"^''"'-'^* - ^ {pq\V''\jr)alalaraka^e" 



pqr 



pqr 



We proceed by dividing the system into two subspaces, namely low (L) and high-energy states (H). The high-energy 
states evolve on a time scale which is short compared to the characteristic evolution at low energy. We can therefore 
adiabatically eliminate operators which act in the high-energy subspace by expressing them in terms of operators which 
act only on low-energy states. This results in the replacement of the bare interaction potential between particles in 
the low-energy subspace with an effective interaction which has the form of a two-body T-matrix. 

In order to prove this result we need to integrate the equation of motion for the product of three operators using 
Eq. (A5) and substitute the res ult into Eq. (A4). Of course this is not possible in general because we do not have 
a closed set of equations. Eq. (A5) depends on the equation of motion for the product of five operators which in 
turn depends on products of seven operators and so on ad infinitum. This infinite set of equations can be truncated 
to just Eqs. (A4) and (A5), however, on the basis of two assumptions. Since we are interested in the properties of 
Bose gases at very low temperatures, we assume that all the particles are confined to the low-energy subspace. The 
high-energy states are only occupied in a virtual sense, as the intermediate states in collisions of low-energy particles. 
The dominant terms in any equation of motion therefore have low indices on all operators. We also assume that 
momentum conservation allows us to neglect all terms involving matrix elements where only one index is high (in a 
trap such matrix elements will be non-zero but should still be negligible) . 



Th ese assumptions, and the fact that we are only interested in the case that z is a low index, mean that Eqs. (A4) 
and (A5) can be reduced to 



L H 



jkm 

d{d]dkd„^ ^ 



j km 
H 



dt 



{km\V\pq)a^^apaqe''' 



{km\V''\pq)a^^apaqe^' 



(A6) 
(A7) 



pq 



pq 



(j Low, k,m High) 



where X^pg means that both the indices p and q are low while means that both are high. In these equations, the 
first term on the right hand side gives the dominant contribution since it contains operators which act on populated 
states. The second term in Eq. (A7) contains operators of the same form as those on the left hand side, so this 
equation has a Lippmann-Schwinger form and its solution introduces a T-matrix. 

We will solve Eq. (A7) using an approach based on Fourier transforms. We define the Fourier transform of djdkdm 
by 



fjkm{uj) = 77- / dt a\aka,n{t) e"^* 

f + 00 

d]dkd^{t) ^ I du) fjkm{i^)e' 



(A8) 
(A9) 



Using these definitions in Eq. (A7) and initially considering only the first term on the right hand side leads to the 
result 
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^ {km\V''\pq)e'' 



+00 



duj 



(AlO) 



We can choose the boundary between the high and low-energy states so that the evolution of a low-lying state has 
no significant frequency components in the high-energy subspace. We can therefore neglect the dependence on a;, ujp 
and u>q in the energy denominator of Eq. (AlO) to obtain 



{km\V'^\pq)c 



a-apQq 



-hiujk + iOk) ' 



(All) 



We now postulate that the full solution to Eq. (A7) has this form but with the matrix element {km\V''\pq) replaced 
with {km\TY{\pq) , where Tjj is an operator to be determined. We therefore have the ansatz 



a]akam = ^ {km\TY{\pq)t 



pq 



-h{ujk + ujk) 



(A12) 



Repeating the procedure which led to Eq. (All) we find that this is a solution of Eq. (A?) if Th satisfies the Lippmann- 
Schwinger equation 



H 



1 



-h{ujp + LUq 



-{pq\Tii. 



(A13) 



This is the usual expression for the two-body T-matrix except for the fact that the intermediate states \pq) in the 
collision are restricted to the high-energy subspace. For this reason we have denoted the T-matrix with the subscript 
H. We expect that Th will be approximately equal to if the density and temperature of the system are such that 
the states rapidly become perturbative (i.e. if the boundary between L and H can be taken at low energy). The 
relationship between Th a-nd is discussed below. 

If we now substitute Eq. (A12) into Eq. (A4) we obtain the restricted T-matrix 7h in the equation of motion for 
the annihilation operator a; 



jkm 



(A14) 



This has the same form as Eq. ( |A4D but with a summation which extends only over low-lying states, and a modified 
two-body interaction which is the restricted T-matrix. This equation of motion therefore corresponds to an effective 
Hamiltonian for the low-energy subspace which is 



L ^ L 



aka„ 



(A15) 



ijkm 



Finally we can use a unitary transformation to the basis orthogonal to the condensate which is used in the main text. 
The effective Hamiltonian can therefore be written in the form of Eq. (H) but with the substitutions ^ — > and 
{ij\V'^\km) {ii\TYi\km) . This is the result quoted in Sec. ^ 

The restricted T-matrix Th is related to the two-body T-matrix T2h of Eq. (|8^) by 



T2b(^) = TH + ^rH|M)'^: 



pq 



(4' 



— ^P{pq\T2^{z). 
£9 ) 



(A16) 
for 



This result can be proved by substituting the definition of Th from Eq. (A13), which leads directly to Eq, 
T2h if we assume that z can be neglected for high-energy states and that +S = "This last assumption 
means that terms involving one high and one low index are neglected, which is justified by momentum conservation 
as mentioned above. 

We can solve Eq. ( A16D in the homogeneous limit if we assume that both Th and can be approximated by a 
contact potential. Writing Th = gS{r) and T2b = Uq S{r) this equation becomes (for k,m,z = 0) 
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^0 

1 - C/oOk' 



where uk is defined by 



ax 



K 



mK 



(A17) 



(A18) 



and the integral has a high-energy cut-ofF at a wave vector K. If the introduction of Th is to be useful we need to 
have g KiUo which means that Uq ax ^ 1- Eq. ( A18) therefore gives 



Ka < 1. 



(A19) 



Thus the requirement that Tjj ~ sets an upper limit to the boundary between high and low-energy states. A 
lower limit is provided by the requirement that the approximations we made in eliminating the high-energy states 
are valid. Specifically we require that the only significant population occurs in the low-energy subspace and that the 
evolution in this subspace occurs on a time scale which is long compared to that of the high-energy states. At zero 
temperature the second requirement is the critical factor, and the time scale for the low-energy evolution is determined 
by noUo = [fik^Y /2m. The condition that the high-energy states evolve faster than this is therefore 



K > fco- 



(A20) 



K can simultaneously satisfy the competing conditions set by Eqs. (A19) and (A20) if no? ^ 1, i.e. if the system is 
a dilute gas. 

At finite temperature the lower limit to K is set by the first of the two requirements given above, which can be 
written as 



K > 1/A, 



dB, 



(A21) 



where A^b is the de Broglie wave length. Thus K must simultaneously satisfy Eqs. (A19) and (A21) which can be 



done if a/ X^b ^ 1- This condition is satisfied for a dilute gas even at temperatures above the critical point. 
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FIG. 1. One of the polarization diagrams which contributes to AAiij{ep) of Eq. (BS 
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